
!VERSION NUMBER:
!  $Id$

!module donner_cloud_model_inter_mod

!#include "donner_cloud_model_interfaces.h"

!end module donner_cloud_model_inter_mod



!#####################################################################

subroutine don_cm_cloud_model_k   &
         (nlev_lsm, nlev_hires, ntr, kou, diag_unit, debug_ijt, Param, &
!++lwh
          Col_diag, Initialized, tb, pb, alpp, cld_press, temp_c, &
!--lwh
          mixing_ratio_c, pfull_c, phalf_c, tracers_c, pcsave, &
          exit_flag_c, wv, rcl, dpf, dpftr, qlw, dfr, flux, pt_kou,  &
          dint, cu, cell_precip, apt, cell_melt, pmelt_lsm, summel, &
          efchr, emfhr, cfracice, etfhr, ncc_kou, tcc, ermesg, error)

!--------------------------------------------------------------------
!
!                ONE-DIMENSIONAL CLOUD MODEL 
!                L. DONNER     NCAR     3 OCT 1984
!
!    subroutine cloud_model receives as input cloud base temperature 
!    (tb), pressure (pb), large-scale model vertical profiles of temper-
!    ature (temp_c), mixing ratio (mixing_ratio_c), pressure (pfull_c and
!    phalf_c), and tracer concentrations (tracers_c) for ensemble member
!    kou and produces as output the in-cloud profiles of temperature 
!    (tcc), vertical velocity (wv), cloud radius (rcl), liquid water 
!    (qlw), condensation rate (dpf), wet-deposition rate (dpftr), 
!    freezing rate (dfr), mass flux
!    (flux), tracer concentrations (xclo), the environmental profiles of 
!    temperature (te), mixing ratio (mre) and tracer concentrations 
!    (xtrae), cloud top pressure (pt_kou), and column integrals of pre-
!    cipitation rate (precip), condensation rate (conint) and freezing 
!    rate (dint). if the current column is a column for which diagnostics
!    is desired (debug_ijt = .true.), then output is written to the diag-
!    nostics file (diag_unit).
!--------------------------------------------------------------------

use donner_types_mod, only : donner_param_type, donner_column_diag_type, &
!++lwh
                             donner_initialized_type
!--lwh

implicit none

!--------------------------------------------------------------------
integer,                            intent(in)     :: nlev_lsm,  &
                                                      nlev_hires, ntr, &
                                                      kou, diag_unit
logical,                            intent(in)     :: debug_ijt
type(donner_param_type),            intent(in)     :: Param
type(donner_column_diag_type),      intent(in)     :: Col_diag
!++lwh
type(donner_initialized_type),      intent(in)     :: Initialized
!--lwh
real,                               intent(in)     :: tb, pb, alpp
real,    dimension(nlev_hires),     intent(in)     :: cld_press
real,    dimension(nlev_lsm),       intent(in)     :: temp_c,   &
                                                      mixing_ratio_c, &
                                                      pfull_c
real,    dimension(nlev_lsm+1),     intent(in)     :: phalf_c 
real,    dimension(nlev_lsm,ntr),   intent(in)     :: tracers_c
real,                               intent(inout)  :: pcsave
logical,                            intent(inout)  :: exit_flag_c     
real,    dimension(nlev_hires),     intent(out)    :: wv, rcl, dpf, &
                                                      qlw, dfr, flux, &
                                                      tcc 
real,                               intent(out)    :: pt_kou, dint, cu,&
                                                      cell_precip,  &
                                                             apt
real,    dimension(nlev_lsm),       intent(out)    :: cell_melt
real,    dimension(nlev_hires),     intent(out)    :: efchr, emfhr,   &
                                                      cfracice
real,    dimension(nlev_hires,ntr), intent(out)    :: etfhr, dpftr
integer,                            intent(out)    :: ncc_kou
real,                               intent( in)    :: pmelt_lsm
real,                               intent(out)    :: summel
character(len=*),                   intent(out)    :: ermesg
integer,                            intent(out)    :: error

!---------------------------------------------------------------------
!   intent(in) variables:
!
!     tb             cloud base temperature [ deg K ]
!     pb             cloud base pressure [ Pa ]             
!     temp_c         large-scale model temperature profile (index 1 
!                    nearest the surface) [ deg K ]
!     mixing_ratio_c large-scale model mixing ratio profile (index 1 
!                    nearest the surface) [ kg(h2o)/ kg(air) ]
!     sig            sigma coordinate of large-scale model levels
!                    index 1 nearest the surface) [ dimensionless ]
!     tracers_c      large-scale model tracer concentration profiles 
!                    (index 1 nearest the surface) [ kg/ kg ]  ???
!     kou            current ensemble member index i
!     diag_unit      output unit number for this diagnostics column
!     debug_ijt      is this a diagnostics column ?
!
!   intent(out) variables:
!
!     tcc            in-cloud temperature profile (index 1 at physical 
!                    base of cloud)  [ deg K ]
!     wv             in-cloud vertical velocity profile (index 1 at 
!                    physical base of cloud)  [ m / sec ]
!     rcl            cloud radius profile (index 1 at physical 
!                    base of cloud)  [ m ]
!     dpf            cloud-area-weighted condensation rate profile
!                    (index 1 at physical base of cloud) 
!                    [ (m**2) * kg(h2o) / kg(air) / sec ]
!     dpftr          wet-deposition rate profile,
!                    weighted by ratio of fractional area to
!                    fractional area at base 
!                    (index 1 at physical base of cloud)
!                    [  [units of xclo] /(sec) ]
!     qlw            cloud liquid water content profile (index 1 at
!                    physical base of cloud) [ kg(h2o) / kg(air) ]
!     dfr            cloud-area-weighted freezing rate profile in 
!                    convective updraft (index 1 at physical base of 
!                    cloud) [ (m**2) *g(h2o) / (kg(air) /  day ]
!     flux           upward mass flux profile in cloud (index 1 at 
!                    physical base of cloud) [ kg(air) / sec ]
!     xclo           in-cloud tracer concentration profiles (index 1 at 
!                    physical base of cloud)  [ kg / kg ] ???
!     te             environmental temperature profile on cloud-model 
!                    grid (index 1 at physical base of cloud) [ deg K ]
!     mre            environmental mixing ratio profile on cloud-model 
!                    grid (index 1 at physical base of cloud) 
!                    [ kg(air) / kg(h2o) ]
!     cell_melt      in-cloud melting of condensate associated with 
!                    convective cells. made up of two parts, 1) that 
!                    due to the freezing of liquid carried upwards 
!                    in the cell updraft, 2) that due to the melting of
!                    condensed ice that precipitates out. if meso 
!                    circulation is present, this component is zero; 
!                    melting will be determined in subroutine mesub.
!     xtrae          environmental tracer profiles on cloud-model 
!                    grid (index 1 at physical base of cloud) 
!                    [ kg / kg ]  ?????
!     pt_kou         pressure at cloud top for this ensemble member[ Pa ]
!     precip         total precipitation rate [ kg(h2o) / sec ]
!     conint         total condensation rate [ kg(h2o) / sec ]
!     dint           total freezing rate in convective updraft 
!                    [ kg(h20) / sec ]
!     cu
!     rc
!     lfc_not_reached      level of free convection was reached ?
!
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!   local variables:

      real, dimension (nlev_hires)         ::  pf, te, mre
      real, dimension (nlev_hires,ntr)     ::  xclo, xtrae, pftr 
      real        :: precip, conint,                 pmel
      integer     :: k, kc
      real        :: accond, acpre
      real        :: sumfrea, sumlhr
      logical     :: lfc_not_reached, do_donner_tracer
      integer     ::  cldtop_indx

!--------------------------------------------------------------------
!   local variables:
!
!
!--------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!--------------------------------------------------------------------
!    if in diagnostic column, output the ensemble member index which
!    is being integrated.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, i4)')    &
            'PROCESSING CLOUD ENSEMBLE MEMBER # ', kou
      endif

!-------------------------------------------------------------------
!    define  or initialize an array to hold any in-cloud tracer sources.
!-------------------------------------------------------------------
      call don_cm_gen_incloud_profs_k  &
           (nlev_lsm, nlev_hires, ntr, kou, diag_unit, debug_ijt, &
            Col_diag, Param, Initialized, tb, pb, alpp, cld_press,  &
            temp_c, mixing_ratio_c, pfull_c, phalf_c, tracers_c,  &
            pcsave,tcc, wv, rcl, qlw, dfr, flux, pf, pftr, te, mre, &
            xclo, xtrae, dint, accond, acpre, sumfrea, cldtop_indx, &
            do_donner_tracer, lfc_not_reached, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!--------------------------------------------------------------------
!    define the cloud top pressure for the current ensemble member.
!--------------------------------------------------------------------
      pt_kou = pb + real(cldtop_indx)*Param%dp_of_cloud_model

!--------------------------------------------------------------------
!    if in diagnostics column, output the ensemble member number and
!    the level of free convection of the most-entraining (first)
!    ensemble member.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, i4, f19.10)') &
                        'in cloudm: kou,pcsave= ',kou,pcsave
      endif

!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
      call don_cm_process_condensate_k     &
         (nlev_lsm, nlev_hires, ntr, cldtop_indx, diag_unit, debug_ijt,&
          Param, acpre, accond, pb, pt_kou, pf, pftr, tcc, rcl,  &
          cld_press, phalf_c, conint, dint, pmel, pmelt_lsm, precip, &
          cu, cell_precip, sumlhr, summel, dpf, dpftr, dfr, cell_melt, &
          ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!----------------------------------------------------------------------
!    if the level of free convection was never reached for this ensemble
!    member, set exit_flag_c so that calculations will cease in this 
!    column, and exit the ensemble members loop. if this is a diag-
!    nostics column, print a message.
!----------------------------------------------------------------------
      if (lfc_not_reached) then
        if (debug_ijt) then
          write (diag_unit, '(a)')  &
                           'in mulsub: lfc never reached'     
        endif
        exit_flag_c = .true.
        return
      endif

!----------------------------------------------------------------------
!    if no precipitation was produced by the ensemble member, set 
!    exit_flag_c so that calculations will cease in this column, and exit
!    the ensemble members loop. if this is a diagnostics column, print 
!    a message.
!----------------------------------------------------------------------
      if (precip == 0.) then
        if (debug_ijt) then
          write (diag_unit, '(a)') &
                       'in mulsub: PRECIP=0 AFTER CLOUD MODEL'
        endif
        exit_flag_c   = .true.
        return
      endif

!---------------------------------------------------------------------
!    if condensate is being evaporated at any level within the cloud, 
!    set exit_flag_c so that calculations will cease in this 
!    column, and exit the ensemble members loop. if this is a diag-
!    nostics column, print a message.
!---------------------------------------------------------------------
      do kc=1,nlev_hires-1           
        if (dpf(kc) > 0.) then 
          if (debug_ijt) then
            write (diag_unit, '(a)') 'in mulsub: dpf  .GT. 0.'
          endif
          exit_flag_c   = .true.
          return   
        endif 
      end do
      if (exit_flag_c) return

!--------------------------------------------------------------------
!    define the cloud model vertical index that is just above cloud top
!    (ncc_kou). this index may be used to limit calculations on the cloud
!    model dimensioned arrays.
!--------------------------------------------------------------------
      do k=1,nlev_hires
        if (cld_press(k) < pt_kou) then
          ncc_kou = k 
          exit
        endif
      end do

!--------------------------------------------------------------------
!    call compute_vertical_fluxes to calculate cumulus thermal forcing 
!    and moisture forcing on the cloud-model grid associated with this
!    ensemble member. 
!--------------------------------------------------------------------
      call don_cm_compute_vert_fluxes_k   &
           (nlev_hires, ntr, ncc_kou, kou, diag_unit, debug_ijt, &
            do_donner_tracer, Param, pt_kou, cld_press, rcl, te, mre, &
            wv, tcc, dpf, xclo, xtrae, apt, efchr, emfhr, etfhr, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!----------------------------------------------------------------------
!    define the fraction of condensation that is ice (cfraci) and that
!    which is liquid (cfracl). at temps above tfre, all condensate is
!    liquid;at temps below (tfre-dfre) all condensate is ice; in between
!    the condensate is partitioned linearly with the temperature depart-
!    ure between these two limits.
!----------------------------------------------------------------------
      cfracice = 0.
      do k=1,ncc_kou
        if (qlw(k) > 0.0) then
          if (tcc(k) > Param%tfre) then
            cfracice(k) = 0.
          else if (tcc(k) < (Param%tfre - Param%dfre)) then
            cfracice(k) = 1.
          else
            cfracice(k) = (Param%tfre - tcc(k))/Param%dfre
          endif
        endif
      end do

!---------------------------------------------------------------------


end subroutine don_cm_cloud_model_k



!#####################################################################

subroutine don_cm_gen_incloud_profs_k  &
         (nlev_lsm, nlev_hires, ntr, kou, diag_unit, debug_ijt, &
          Col_diag,  Param, Initialized, tb, pb, alpp, cld_press, &
          temp_c, mixing_ratio_c, pfull_c, phalf_c, tracers_c, pcsave, &
          tcc, wv, rcl, qlwa, dfr, flux, pf, pftr, te, mre, xclo, &
          xtrae, dint, accond, acpre, sumfrea, cldtop_indx,  &
          do_donner_tracer, lfc_not_reached, ermesg, error)

!----------------------------------------------------------------------
!
!----------------------------------------------------------------------

!  modified by Leo Donner, GFDL, 5 February 2007
!
use donner_types_mod, only : donner_param_type, &
                             donner_column_diag_type, &
                             donner_initialized_type
use sat_vapor_pres_k_mod, only: compute_mrs_k

implicit none 

!----------------------------------------------------------------------
integer,                            intent(in)    :: nlev_lsm,   &
                                                     nlev_hires,&
                                                     ntr, kou, diag_unit
logical,                            intent(in)    :: debug_ijt
type(donner_column_diag_type),      intent(in)    :: Col_diag
type(donner_param_type),            intent(in)    :: Param      
!++lwh
type(donner_initialized_type),      intent(in)    :: Initialized
!--lwh
real,                               intent(in)    :: tb, pb, alpp
real,    dimension(nlev_hires),     intent(in)    :: cld_press
real,    dimension(nlev_lsm),       intent(in)    :: temp_c,   &
                                                     mixing_ratio_c,  &
                                                     pfull_c
real,    dimension(nlev_lsm+1),     intent(in)    :: phalf_c    
real,    dimension(nlev_lsm,ntr),   intent(in)    :: tracers_c
real,                               intent(inout) :: pcsave         
real,    dimension(nlev_hires),     intent(out)   :: tcc, wv, rcl, qlwa,&
                                                     dfr, flux, pf
real,    dimension(nlev_hires),     intent(out)   :: te, mre
real,    dimension(nlev_hires,ntr), intent(out)   :: xclo, xtrae, pftr
real,                               intent(out)   :: dint, accond, acpre
real,                               intent(out)  :: sumfrea
integer,                            intent(out)   :: cldtop_indx
logical,                            intent(out)   :: do_donner_tracer, &
                                                     lfc_not_reached
character(len=*),                   intent(out)   :: ermesg
integer,                            intent(out)   :: error

!--------------------------------------------------------------------
!
!--------------------------------------------------------------------


!---------------------------------------------------------------------
!   intent(in) variables:
!
!     psfc           surface pressure [ Pa ]
!     tb             cloud base temperature [ deg K ]
!     pb             cloud base pressure [ Pa ]             
!     temp_c         large-scale model temperature profile (index 1 
!                    nearest the surface) [ deg K ]
!     mixing_ratio_c large-scale model mixing ratio profile (index 1 
!                    nearest the surface) [ kg(h2o)/ kg(air) ]
!     sig            sigma coordinate of large-scale model levels
!                    index 1 nearest the surface) [ dimensionless ]
!     tracers_c      large-scale model tracer concentration profiles 
!                    (index 1 nearest the surface) [ kg/ kg ]  ???
!     kou            current ensemble member index i
!     diag_unit      output unit number for this diagnostics column
!     debug_ijt      is this a diagnostics column ?
!
!   intent(out) variables:
!
!     tcc            in-cloud temperature profile (index 1 at physical 
!                    base of cloud)  [ deg K ]
!     wv             in-cloud vertical velocity profile (index 1 at 
!                    physical base of cloud)  [ m / sec ]
!     rcl            cloud radius profile (index 1 at physical 
!                    base of cloud)  [ m ]
!     pf             cloud-area-weighted condensation rate profile
!                    (index 1 at physical base of cloud) 
!                    [ (m**2) * kg(h2o) / kg(air) / sec ]
!     pftr           cloud-area-weighted wet-deposition profile
!                    (index 1 at physical base of cloud)
!                    [ (m**2) * [xlco units] /sec ]
!     qlwa           cloud liquid water content profile (index 1 at
!                    physical base of cloud) [ kg(h2o) / kg(air) ]
!     dfr            cloud-area-weighted freezing rate profile in 
!                    convective updraft (index 1 at physical base of 
!                    cloud) [ (m**2) *g(h2o) / (kg(air) /  day ]
!     flux           upward mass flux profile in cloud (index 1 at 
!                    physical base of cloud) [ kg(air) / sec ]
!     xclo           in-cloud tracer concentration profiles (index 1 at 
!                    physical base of cloud)  [ kg / kg ] ???
!     te             environmental temperature profile on cloud-model 
!                    grid (index 1 at physical base of cloud) [ deg K ]
!     mre            environmental mixing ratio profile on cloud-model 
!                    grid (index 1 at physical base of cloud) 
!                    [ kg(air) / kg(h2o) ]
!     cell_melt      in-cloud melting of condensate associated with 
!                    convective cells. made up of two parts, 1) that 
!                    due to the freezing of liquid carried upwards 
!                    in the cell updraft, 2) that due to the melting of
!                    condensed ice that precipitates out. if meso 
!                    circulation is present, this component is zero; 
!                    melting will be determined in subroutine mesub.
!     xtrae          environmental tracer profiles on cloud-model 
!                    grid (index 1 at physical base of cloud) 
!                    [ kg / kg ]  ?????
!     pt             pressure at cloud top [ Pa ]
!     precip         total precipitation rate [ kg(h2o) / sec ]
!     conint         total condensation rate [ kg(h2o) / sec ]
!     dint           total freezing rate in convective updraft 
!                    [ kg(h20) / sec ]
!     rc
!     lfc_not_reached      level of free convection was reached ?
!
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!   local variables:

      real,   dimension(nlev_lsm)         :: sig
      real,   dimension(nlev_hires)       :: dis, rsc
      real,   dimension (nlev_hires,ntr)  ::   clsou
      real     :: qlw, qcw, qrw, dtfr, dtupa, dfrac, rmu,  &
                  rhodt_inv, psfc, actot, dt_micro,          dcw1,  &
                  dqrw3, rbar, rmub, density, densityp, d1, d2, dztr, &
                  entrain, dt_inv
      logical  :: flag
      integer  :: max_cloud_level, ktr
      integer  :: k, nbad
!++lwh
      real, parameter :: g_2_kg = 1.e-3 ! kg/g
      real :: qlw_save, t_avg
      real, dimension( size(xclo,2) ) :: delta_xclo0, delta_xclo1, dwet
      integer :: n
!--lwh
      real :: displ
      integer :: indx

!--------------------------------------------------------------------
!   local variables:
!
!
!--------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = '  ' ; error = 0

      psfc = phalf_c(1)
      if (ntr /= 0) then
        do_donner_tracer = .true.
      else
        do_donner_tracer = .false.
      endif

!--------------------------------------------------------------------
!    provide appropriate initialization for the output fields. 
!    NOTE: input tracer fields of all zeroes are currently provided if 
!    do_donner_tracer is .false.. one could make the tracer fields
!    optional and define output only when present.
!--------------------------------------------------------------------
      tcc(:)  = 0.
      wv(:)   = 0.
      rcl(:)  = 0.
      qlwa(:) = 0.
      dfr(:)  = 0.
      te(:)   = 0.
      mre(:)  = 0.
      xclo(:,:)  = 0.
      xtrae(:,:) = 0.
      pftr(:,:) = 0.
      flux(:) = 0.
      dint   = 0. 
      lfc_not_reached = .true.

!-------------------------------------------------------------------
!    initialize local variables.
!-------------------------------------------------------------------
      do k=1,nlev_hires
        dis(k)=0.
        pf(k)=0.
      end do

      rsc(:) = 0.
      sig(:) = pfull_c(:)/psfc

!--------------------------------------------------------------------
!    if in diagnostics column, output the cloud base radius, pressure
!    and temperature, and the large-scale model sigma levels.
!--------------------------------------------------------------------
      if (debug_ijt) then
        do k=1,nlev_lsm-Col_diag%kstart+1
          write (diag_unit, '(a, i4, e20.12)')  &
                        'in cloudm: k,sig   = ',k, sig(k)    
        end do
      endif

!---------------------------------------------------------------------
!    initialize various integrals.
!---------------------------------------------------------------------

      qcw = 0.
      qlw = 0.
      qrw = 0.

!---------------------------------------------------------------------
!    initialize various integrals.
!---------------------------------------------------------------------
      accond  = 0.
      acpre   = 0.
      dtupa   = 0.
      dfrac   = 0.
      dtfr    = 0.
      sumfrea = 0.

!---------------------------------------------------------------------
!    initialize various conditionals.
!---------------------------------------------------------------------
      cldtop_indx = 0
  
!--------------------------------------------------------------------
!    define the cloud model pressure levels. 
!--------------------------------------------------------------------
      max_cloud_level = nlev_hires + 2
      do k=1,nlev_hires
        if (cld_press(k) < Param%pstop) then
          max_cloud_level = k - 1 
          exit
        endif
      end do

!--------------------------------------------------------------------
!    define cloud base (k = 1) values of vertical velocity (wv),
!    temperature (tcc), mixing ratio (rsc), cloud radius (rcl) and
!    tracer concentrations (xclo).
!--------------------------------------------------------------------
      wv(1) = Param%cld_base_vert_vel
      rcl(1) = Param%cloud_base_radius
      tcc(1) = tb
      xclo(1,:) = tracers_c(1,:)
      call compute_mrs_k (tb, pb, Param%D622, Param%D608, rsc(1), nbad)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (nbad /= 0) then
        ermesg = 'subroutine don_cm_gen_incloud_profs_k: '// &
                 'temperatures out of range of esat table'
        error = 1
        return
      endif

!--------------------------------------------------------------------
!    call gcm_to_cm to obtain the large-scale model moisture, temper-
!    ature, and tracer field values at the cloud base pressure to be
!    used as the environmental values in the cloud model.
!
!    call don_u_lo1d_to_hi0d_log_fast_k to determine the index and 
!    displacement needed to be used on all interpolated variables 
!    (mre, te, tracers)
!--------------------------------------------------------------------
      call don_u_lo1d_to_hi0d_log_fast_k             &
           (nlev_lsm, sig, psfc, pb,  &
            indx, displ, ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

      mre(1) = mixing_ratio_c(indx) + &
              (mixing_ratio_c(indx+1) - mixing_ratio_c(indx))*displ
      te(1) = temp_c(indx) + &
              (temp_c(indx+1) - temp_c(indx))*displ

      if (do_donner_tracer) then
        do ktr=1,ntr          
          xtrae(1,ktr) = tracers_c(indx,ktr) + &
                (tracers_c(indx+1,ktr) - tracers_c(indx,ktr))*displ
        end do
      endif

!--------------------------------------------------------------------
!    if in diagnostics column, output the cloud base radius, pressure
!    and temperature, and the large-scale model sigma levels.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12, f19.10, f20.14)')  &
                'in cloudm: RR,PB,TB= ',Param%cloud_base_radius, pb, tb
      endif

!--------------------------------------------------------------------
!    if in diagnostics column, output cloud base values of environ-
!    mental moisture (mre) and temperature (te), and the in-cloud 
!    temperature (tcc).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12,f20.14)')  &
                       'in cloudm: QE,TE,TCC= ',mre(1), te(1), tcc(1)
      endif
      
!---------------------------------------------------------------------
!    loop over cloud model levels, generating the desired in-cloud and 
!    environmental profiles. under no circumstances extend the cal-
!    culation above max_cloud_level (the level above which cloud is not
!    allowed).
!---------------------------------------------------------------------
      do k=1,max_cloud_level-1

!----------------------------------------------------------------------
!    exit the loop if the pressure is lower than the upper limit allowed
!    for the level of free convection and free convection has not yet 
!    been achieved.
!----------------------------------------------------------------------
        if ((cld_press(k+1) <= Param%upper_limit_for_lfc) .and. &
            (lfc_not_reached)) then
          cldtop_indx = 0
          exit
        endif 

!--------------------------------------------------------------------
!    call gcm_to_cm to obtain the large-scale model moisture, temper-
!    ature, and tracer field values at cloud model pressure level to be
!    used as the environmental values in the cloud model.
!
!    call don_u_lo1d_to_hi0d_log_fast_k to determine the index and 
!    displacement needed to be used on all interpolated variables 
!    (mre, te, tracers)
!--------------------------------------------------------------------
        call don_u_lo1d_to_hi0d_log_fast_k             &
             (nlev_lsm, sig, psfc, cld_press(k+1),  &
              indx, displ, ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (error /= 0 ) return

        mre(k+1) = mixing_ratio_c(indx) + &
                  (mixing_ratio_c(indx+1) - mixing_ratio_c(indx))*displ

        te(k+1) = temp_c(indx) + &
                 (temp_c(indx+1) - temp_c(indx))*displ

        if (do_donner_tracer) then
          do ktr=1, ntr           
            xtrae(k+1,ktr) = tracers_c(indx,ktr) + &
                            (tracers_c(indx+1,ktr) - tracers_c(indx,ktr))*displ
          end do
        endif

!--------------------------------------------------------------------
!    if in diagnostics column, output level k values of vertical
!    velocity, pressure, environmental temperature (te) and moisture 
!    (mre), the in-cloud temperature (tcc) and liquid water (qlw).
!--------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, 3e20.12)')  &
                      'in cloudm: WV,PP,TE= ',wv(k), cld_press(k), te(k)
          write (diag_unit, '(a,f20.14, 2e20.12)')  &
                     'in cloudm: TCC(k),qe(k),QLW= ',tcc(k), wv(k), qlw
        endif

!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!++lwh
        qlw_save = qlw
!--lwh

        call don_cm_move_parcel_k    &
             (k, kou, diag_unit, debug_ijt, cld_press(k),  &
              cld_press(k+1), alpp, Param, pcsave, qlw, sumfrea, qrw, &
              qcw, dcw1, dqrw3, wv(k), tcc(k), rcl(k), te(k), mre(k), &
              rcl(k+1), wv(k+1), tcc(k+1), rsc(k+1), te(k+1), mre(k+1),&
              qlwa(k+1), dfr(k+1), rmu, rbar, dfrac, dtfr, dtupa,  &
              lfc_not_reached, flag, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (error /= 0 ) return

        if (flag) exit

!---------------------------------------------------------------------
!    define the in-cloud  density at levels k and k+1 
!    (density, densityp).
!---------------------------------------------------------------------
        density  = cld_press(k)/(Param%rdgas*(tcc(k)*  &
                               (1. + Param%D608*(rsc(k)/(1.0+rsc(k))))))
        densityp = cld_press(k+1)/(Param%rdgas*(tcc(k+1)*  &
                           (1. + Param%D608*(rsc(k+1)/(1.0+rsc(k+1))))))
 
!--------------------------------------------------------------------
!    calculate in-cloud tracer distribution.
!--------------------------------------------------------------------
        if (do_donner_tracer) then
          clsou(k,:) = 0.0
          clsou(k+1,:) = 0.0
          d1 = Param%rdgas*(1.+Param%virt_mass_co)*tcc(k)*te(k)/  &
               (Param%grav*cld_press(k)*(Param%virt_mass_co*te(k) +  &
                                                                tcc(k)))
          d2 = Param%rdgas*(1.+Param%virt_mass_co)*tcc(k+1)*te(k+1)/ &
               (Param%grav*cld_press(k+1)*      &
                                (Param%virt_mass_co*te(k+1) + tcc(k+1)))
          dztr = ((d1 + d2)*(cld_press(k) - cld_press(k+1))/2.)
          rmub = 2.*alpp/rbar
          entrain = dztr*rmub
          dt_micro = dztr/(0.5*(wv(k) + wv(k+1)))
          if (debug_ijt) then
            write (diag_unit, '(a, 3e20.12)')   &
                 'in clotr: d1,d2,dz= ',d1,d2,dztr
            write (diag_unit, '(a, e20.12)')   &
                 'in clotr: ent= ',rmub
          endif

!++lwh
!--------------------------------------------------------------------
!    Call tracer wet deposition
!    
!    Convert dqrw3 from g/m3 to kg/m3
!    from g(h2o) per m**3 to kg(h2o) per kg(air).
!--------------------------------------------------------------------
          t_avg = 0.5*(tcc(k)+tcc(k+1))
          do n = 1,size(xclo,2)
             if (Initialized%wetdep(n)%Lwetdep) then
                call wet_deposition_0D    &
                           (Initialized%wetdep(n)%Henry_constant, &
                            Initialized%wetdep(n)%Henry_variable, &
                            Initialized%wetdep(n)%frac_in_cloud, &
                            Initialized%wetdep(n)%alpha_r, &
                            Initialized%wetdep(n)%alpha_s, &
                            t_avg, cld_press(k), cld_press(k+1), &
                            0.5*(density+densityp), &
                            qlw_save, dqrw3*g_2_kg, 0., &
                            xclo(k,n), &
                            Initialized%wetdep(n)%Lgas, &
                            Initialized%wetdep(n)%Laerosol, &
                            Initialized%wetdep(n)%Lice, &
                            delta_xclo0(n) )
             end if
          end do
!--lwh          
          call don_cm_clotr_k    &
               (ntr, diag_unit, debug_ijt, Param, clsou(k,:),  &
                clsou(k+1,:), xtrae(k,:), xtrae(k+1,:), xclo(k,:), &
                entrain, dt_micro, xclo(k+1,:), ermesg, error)
!++lwh
          do n = 1,size(xclo,2)
             if (Initialized%wetdep(n)%Lwetdep) then
                call wet_deposition_0D   &
                           (Initialized%wetdep(n)%Henry_constant, &
                            Initialized%wetdep(n)%Henry_variable, &
                            Initialized%wetdep(n)%frac_in_cloud, &
                            Initialized%wetdep(n)%alpha_r, &
                            Initialized%wetdep(n)%alpha_s, &
                            t_avg, cld_press(k), cld_press(k+1), &
                            0.5*(density+densityp), &
                            qlw, dqrw3*g_2_kg, 0., &
                            xclo(k+1,n), &
                            Initialized%wetdep(n)%Lgas, &
                            Initialized%wetdep(n)%Laerosol, &
                            Initialized%wetdep(n)%Lice, &
                            delta_xclo1(n) )
                dwet(n) = - 0.5*(delta_xclo0(n)+delta_xclo1(n))
                xclo(k+1,n) = xclo(k+1,n) + dwet(n)
             else
                dwet(n) = 0.
             end if
          end do
!--lwh          
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
          if (error /= 0 ) return
        endif

!--------------------------------------------------------------------
!    define cloud-area-weighted precip removal from the layer (dis) in 
!    units of  m**2 * kg(h2o) per kg(air) per day. define rhodt_inv 
!    (10e-03*g*w/delta p) to be used to convert units of dqrw3 and dcw1
!    from g(h2o) per m**3 to kg(h2o) per kg(air).
!--------------------------------------------------------------------
        rhodt_inv = 1.0e-03*(0.5*(wv(k) + wv(k+1)))*Param%GRAV  /   &
                    Param%dp_of_cloud_model
        dis(k) = dqrw3*(rbar**2)*rhodt_inv
        pf(k)  = dcw1*(rbar**2)*rhodt_inv

!---------------------------------------------------------------------
!     define the removal by wet deposition between levels k and k+1
!---------------------------------------------------------------------
        dt_inv = -.25*(wv(k)+wv(k+1))*(density+densityp)*Param%grav  / &
                 Param%dp_of_cloud_model
        do n = 1,size(xclo,2)
          pftr(k,n) = dwet(n) *(rbar**2)*dt_inv
        end do

!---------------------------------------------------------------------
!    calculate moisture subject to freezing in units of g(h2o) per 
!    kg(air) per day, weighted by the cloud area (dfr). add this layer's
!    contribution to the integrated water mass frozen in the updraft 
!    in units of kg(h2o)/sec (dint). 
!---------------------------------------------------------------------
        if (dfr(k+1) /= 0.) then
          dfr(k+1) = dfr(k+1)*densityp*wv(k+1)*(rcl(k+1)**2)*   &
                                      Param%grav/Param%dp_of_cloud_model 
          dfr(k+1) = -dfr(k+1)*Param%cp_air/Param%hlf   
          dint = dint - dfr(k+1)*Param%dp_of_cloud_model/Param%grav  
          dfr(k+1) = dfr(k+1)*8.64E07
        endif 
 
!---------------------------------------------------------------------
!    add this layer's contribution to the integrals accumulating total
!    condensation (accond), precipitation (acpre) and  total liquid
!    water (actot).
!---------------------------------------------------------------------
        accond = accond + dcw1
        acpre  = acpre  + dqrw3
        actot =  accond + acpre

!--------------------------------------------------------------------
!    if in diagnostics column, output the values of accumulated conden-
!    sation, precipitation and total liquid water at this level.
!--------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, i4, 2e20.12)')  &
                        'in cloudm: k,ACCOND,ACPRE= ',k, accond, acpre
          write (diag_unit, '(a, i4, e20.12)')  &
                        'in cloudm:  k,ACTOT= ',k, actot
        endif

!--------------------------------------------------------------------
!    define the mass flux at level k+1 (flux).
!--------------------------------------------------------------------
        if (k == 1) flux(k) = (rcl(k)**2)*wv(k)*density 
        flux(k+1) = (rcl(k+1)**2)*wv(k+1)*densityp

!--------------------------------------------------------------------
!    if in diagnostics column, output the values of cloud radius and
!    entrainment coefficient at this level.
!--------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, i4, 2e20.12)')  &
                     'in cloudm: k,rcl,RMU= ',k, rcl(k), rmu
        endif

!---------------------------------------------------------------------
!    if non-negative water is present, save the current k index as the
!    index of the last cloud model level containing cloud (cldtop_indx).
!---------------------------------------------------------------------
        if (qlwa(k+1) >= 0.) then    
          cldtop_indx = k
        endif  ! (qlwa(k+1) >= 0.0)

!---------------------------------------------------------------------
!    end of loop over cloud model levels. upon leaving loop, in-cloud
!    profiles of model variables have been generated.
!---------------------------------------------------------------------
      end do   ! (do 1 loop)


!--------------------------------------------------------------------
!    if in diagnostics column, output the water mass frozen in the
!    updraft in units of cloud area*kg(h2o) / sec (dint) and the 
!    resultant temperature change (sumfrea).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)')  &
                'in cloudm: DINT IN CLOUDM,sumfrea= ', dint, sumfrea
      endif

!--------------------------------------------------------------------
!    normalize the column integral of water mass frozen in the updraft
!    by the cloud base area of ensemble member #1.
!--------------------------------------------------------------------
      dint = dint/(Param%cloud_base_radius**2)

!--------------------------------------------------------------------
!    if an acceptable cloud was found for this ensemble member, obtain
!    the needed diagnostic / integral output.
!--------------------------------------------------------------------
      if (cldtop_indx /= 0) then

!--------------------------------------------------------------------
!    if in diagnostics column, output the topmost cloud level and
!    the mixing ratios at and one level above cloud top.  
!--------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, i4, 2e20.12)')  &
                      'in cloudm: n,rsc(n+2),rsc(n)= ',cldtop_indx,   &
                             rsc(cldtop_indx+1),rsc(cldtop_indx)
        endif
      endif ! (cldtop_indx /= 0)

!--------------------------------------------------------------------
!    if in diagnostics column, output the values of cloud-area-
!    weighted layer-mean condensation and  fallout.
!--------------------------------------------------------------------
        do k = 1, cldtop_indx
          if (debug_ijt) then
            write (diag_unit, '(a, i4, 2e20.12)')   &
                       'in cloudm: K,PF,DIS= ', k, pf(k), dis(k)
          endif
        end do


end subroutine don_cm_gen_incloud_profs_k


!#####################################################################

subroutine don_cm_move_parcel_k    &
              (k, kou, diag_unit, debug_ijt, pbot, ptop, alpp, Param,  &
               pcsave, qlw, sumfrea, qrw, qcw, dcw1, dqrw3, wvbot,  &
               tccbot, rclbot, tebot, mrebot, rcltop, wvtop, tcctop, &
               rsctop, tetop, mretop, qlwatop, dfrtop, rmu, rbar, &
               dfrac, dtfr,  dtupa, lfc_not_reached, flag, ermesg, error)

!----------------------------------------------------------------------
!
!----------------------------------------------------------------------

use donner_types_mod, only : donner_param_type
use sat_vapor_pres_k_mod, only: compute_mrs_k

implicit none

!----------------------------------------------------------------------
integer,                 intent(in)    :: k, kou, diag_unit
logical,                 intent(in)    :: debug_ijt
real,                    intent(in)    :: pbot, ptop, alpp
type(donner_param_type), intent(in)    :: Param
real,                    intent(inout) :: pcsave, qlw, sumfrea, qrw, &
                                          qcw, dcw1, dqrw3, wvbot,  &
                                          tccbot, rclbot, tebot, &
                                          mrebot, rcltop, wvtop, &
                                          tcctop, rsctop,&
                                          tetop, mretop, qlwatop,  &
                                          dfrtop, rmu, rbar, dfrac,  &
                                          dtfr, dtupa
logical,                 intent(inout) :: lfc_not_reached
logical,                 intent(out)   :: flag
character(len=*),        intent(out)   :: ermesg
integer,                 intent(out)   :: error

!----------------------------------------------------------------------
      real    ::  dtdp, drdp, dwdp, tcest, west, &
                  rest, qcwest,  qlwest, dtdp2, drdp2, dwdp2, qrwest
      integer :: nbad

!---------------------------------------------------------------------
      flag = .false.

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!---------------------------------------------------------------------
!    define the entrainment coefficient (rmu). 
!    pcsave    = p at which ensemble member 1 becomes buoyant.
!    lfc_not_reached = T below level where vert vel starts to increase 
!    with height 
!    THUS:
!    rmu = 0 when above buoyancy level for member 1, above level where 
!            vert vel starts to increase for this member, AND vert vel
!            is so small that are detraining. Thus rmu = 0. when
!            detraining after having been through a convective tower.
!---------------------------------------------------------------------
      if ((pbot    <= pcsave) .and. (.not. lfc_not_reached) .and. &
          (wvbot <= Param%wdet)) then
        rmu = 0.
      else
        rmu = 2.*alpp/rclbot
      end if 

!--------------------------------------------------------------------
!    call simult to solve for the in-cloud derivatives dT/dp (dtdp), 
!    dw/dp (dwdp) and dr/dp (drdp) for a parcel in motion from cloud 
!    model level k (pressure pp) to pressure level pp + delta p.
!--------------------------------------------------------------------
      call don_cm_simult_k   &
           (diag_unit, debug_ijt, lfc_not_reached, Param, pcsave, rmu, &
            tccbot, rclbot, wvbot, pbot, qlw, tebot, mrebot,   &
            dtdp, drdp, dwdp, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!--------------------------------------------------------------------
!    if in diagnostics column, output cloud base values of environ-
!    mental moisture (mre), liquid water (qlw) and initial cloud
!    radius (cloud_base_radius), and the cloud radius and vertical vel-
!    ocity tendencies returned from subroutine simult.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 3e20.12)')   &
             'in cloudm: QE,QLW,RR= ',mrebot,QLW,Param%cloud_base_radius
        write (diag_unit, '(a, 2e20.12)')  &
                                'in cloudm: DPD,DWDP= ',drdp, dwdp
      endif

!--------------------------------------------------------------------
!    define estimated values for t, w and cloud radius based on the
!    derivatives returned from subroutine simult. define the average of
!    the original cloud radius and the new estimate (rbar).
!--------------------------------------------------------------------
      tcest = tccbot + dtdp*Param%dp_of_cloud_model
      west  = wvbot  + dwdp*Param%dp_of_cloud_model
      rest  = rclbot + drdp*Param%dp_of_cloud_model
      rbar  = 0.5*(rest + rclbot)

!--------------------------------------------------------------------
!    if in diagnostics column, output the newly estimated values of 
!    cloud radius (rest) and vertical velocity (west).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)')  &
                         'in cloudm: rest,west= ',rest, west
      endif

!--------------------------------------------------------------------
!    if any of these estimated values are incompatible with the exist-
!    ence of deep convection at this level, exit this vertical loop --
!    cloud does not extend any higher.
!--------------------------------------------------------------------
      if (rbar <= Param%rbound .or. &    
          west <= Param%wbound .or. &     
          rest <= Param%rbound) then
        flag = .true.
        return
      endif

!---------------------------------------------------------------------
!    define estimated values of the liquid water components to be passed
!    to subroutine micro to hold first pass values.
!---------------------------------------------------------------------
      qrwest = qrw
      qcwest = qcw
      qlwest = qlw

!---------------------------------------------------------------------
!    call micro to calculate the microphysical terms which occur during
!    the parcel movement from cld_press(k) to cld_press(k+1).
!---------------------------------------------------------------------
      call don_cm_micro_k   &
           (diag_unit, debug_ijt, Param, tccbot, tcest, pbot, ptop,  &
            tebot, tetop, mrebot, mretop, wvbot, west, rbar, rmu,  &
            qrwest, qcwest, qlwest, dcw1, dqrw3, ermesg, error)  
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!--------------------------------------------------------------------
!    if in diagnostics column, output the newly estimated values of 
!    environmental pressure (cld_press(k+1)), temperature (te) and 
!    mixing ratio (mre), and cloud temperature (tcest), cloud radius 
!    (rest) and liquid water (qlwest)
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, f19.10, f20.14, e20.12)') &
             'in cloudm: P,TEST,QEST= ',ptop     ,TEtop  ,mretop  
        write (diag_unit, '(a, f20.14, 2e20.12)')   &
                   'in cloudm: TCEST,rest,QLWT= ',TCEST,REST,qlwest
      endif

!----------------------------------------------------------------------
!    define the entrainment coefficient (rmu). it is set to 0.0 if the
!    parcel is above the lcl of ensemble member #1, and above the cur-
!    rent ensemble member's cloud base, and if the vertical velocity
!    is below the detrainment threshold; i.e., the current cloud is in
!    its detraining region. if it is not in its detraining zone, set the
!    entrainment coefficient to be the specified entrainment coefficient
!    divided by current cloud radius.
!----------------------------------------------------------------------
      if ((ptop <= pcsave) .and. (.not. lfc_not_reached) .and. &
          (west <= Param%wdet)) then
        rmu = 0.
      else
        rmu = 2.*alpp/rest
      endif 

!----------------------------------------------------------------------
!    call simult to solve for the in-cloud derivatives dT/dp (dtdp2), 
!    dw/dp (dwdp2) and dr/dp (drdp2) using the previously estimated
!    values at p + delta p, including the microphysical contributions.
!----------------------------------------------------------------------
      call don_cm_simult_k   &
           (diag_unit, debug_ijt, lfc_not_reached, Param, pcsave, rmu, &
            tcest, rest, west, ptop, qlwest, tetop, mretop, dtdp2,  &
            drdp2, dwdp2, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!--------------------------------------------------------------------
!    if in diagnostics column, output the lfc flag (lfc_not_reached), and
!    the temperature, cloud radius and vertical velocity tendencies 
!    returned from subroutine simult.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, l4, 2e20.12)')  &
            'in cloudm: TESTLC,DTDP2,DPD2= ',lfc_not_reached,DTDP2,drdp2
        write (diag_unit, '(a, e20.12)')  'in cloudm: DWDP2= ',DWDP2
      endif

!--------------------------------------------------------------------
!    define new values of cloud radius and vertical velocity at level
!    k+1 using the values at level k and the vertical derivatives
!    returned from subroutine simult. if either of these values is below
!    its acceptable bound, set both values to 0.0 and exit the loop -- 
!    the cloud top has been reached.
!--------------------------------------------------------------------
      rcltop   = rclbot + drdp2*Param%dp_of_cloud_model
      wvtop    = wvbot  + dwdp2*Param%dp_of_cloud_model
      if ((wvtop <= Param%wbound) .or. (rcltop   <= Param%rbound)) then 
        rcltop   = 0.0
        wvtop    = 0.
        flag = .true.
        return
      endif

!--------------------------------------------------------------------
!    average the two derivatives calculated for temperature, vertical
!    velocity and cloud radius. define the cloud temperature at the
!    k+1 level as the value at k + the  effect of the calculated deriv-
!    ative.
!--------------------------------------------------------------------
      dtdp = 0.5*(dtdp + dtdp2)
      dwdp = 0.5*(dwdp2 + dwdp)
      drdp = 0.5*(drdp + drdp2)
      tcctop   = tccbot + dtdp*Param%dp_of_cloud_model

!--------------------------------------------------------------------
!    call subroutine freeze_liquid to define the amount of liquid in 
!    the updraft is frozen in the current layer.
!--------------------------------------------------------------------
      call don_cm_freeze_liquid_k   &
           (k, diag_unit, debug_ijt, Param, tccbot, tcctop, qlwest, &
            dfrac, dtfr, dtupa, dfrtop, sumfrea, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!--------------------------------------------------------------------
!    add the effect of droplet freezing to the cloud temperature.
!--------------------------------------------------------------------
      tcctop   = tcctop   + dfrtop  

!--------------------------------------------------------------------
!    define the values of vertical velocity and cloud radius at level 
!    k+1 using the averaged values of dw/dp and dr/dp. if either of 
!    these values is below its acceptable bound, set both values to 0.0
!    and exit the loop -- the cloud top has been reached.
!--------------------------------------------------------------------
      wvtop    = wvbot  + dwdp*Param%dp_of_cloud_model
      rcltop   = rclbot + drdp*Param%dp_of_cloud_model
      if ((wvtop <= Param%wbound) .or. (rcltop <= Param%rbound)) then 
        rcltop   = 0.0
        wvtop    = 0.
        flag = .true.
        return
      endif
 
!--------------------------------------------------------------------
!    define the pressure at which nsemble member #l becomes buoyant (its
!    level of free convection). all ensemble members must have their 
!    cloud top pressure lower than this value.
!--------------------------------------------------------------------
      if (kou == 1) then
        if ((lfc_not_reached) .and. (wvtop > wvbot))  then
          pcsave = ptop     
        endif
      endif

!--------------------------------------------------------------------
!    check to see if the level of free convection has been reached. if
!    so, set the logical variable indicating such (lfc_not_reached) to  
!    be .false..
!--------------------------------------------------------------------
      if (wvtop   > wvbot)  then
        lfc_not_reached = .false.
      endif

!--------------------------------------------------------------------
!    define the vapor mixing ratio for the new cloud temperature at
!    level k+1.
!--------------------------------------------------------------------
      call compute_mrs_k (tcctop, ptop, Param%D622, Param%D608,  &
                          rsctop, nbad)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (nbad /= 0) then
        ermesg = 'subroutine don_cm_move_parcel_k: '// &
                 'temperatures out of range of esat table'
        error = 1
        return
      endif

!--------------------------------------------------------------------
!    if in diagnostics column, output the values of environmental 
!    temperature (te) and mixing ratio (mre) and in-cloud mixing ratio 
!    (rsc) at level k+1.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, f20.14, 2e20.12)')  &
                   'in cloudm: TE,QE,RSC= ',tetop  , mretop  , rsctop  
      endif

!--------------------------------------------------------------------
!    define the layer-mean cloud radius (rbar) and entrainment coef-
!    ficient (rmub).  if the cloud radius is less than the lower bound,
!    exit the loop -- the cloud top has been reached.
!--------------------------------------------------------------------
      rbar = 0.5*(rclbot + rcltop  )
      if (rbar <= Param%rbound) then
        flag = .true.
        return   
      endif

!---------------------------------------------------------------------
!    call micro to calculate the microphysical terms which occur during
!    the parcel movement from p to p + delta p using the second itera-
!    tion values for the cloud conditions at level k + 1.
!---------------------------------------------------------------------
      call don_cm_micro_k   &
           (diag_unit, debug_ijt, Param, tccbot, tcctop, pbot, ptop, &
            tebot, tetop, mrebot, mretop, wvbot, wvtop, rbar, rmu,  &
            qrw, qcw, qlw, dcw1, dqrw3, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!---------------------------------------------------------------------
!    save the final values of cloud water, rainwater and total liquid
!    that are found at level k+1 in arrays disc, disd and qlwa.
!---------------------------------------------------------------------
      qlwatop = qlw

!--------------------------------------------------------------------
!    if in diagnostics column, output the values of in-cloud temperature
!    (tcc), vertical velocity (wv), liquid water (qlw), condensation
!    (dcw1) and precipitation (dqrw3) at level k+1.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, f20.14, 2e20.12)') &
                     'in cloudm: TCC(k+1),WV(k+1),QLW= ',tcctop  ,  &
                      wvtop  , qlw
        write (diag_unit, '(a, 2e20.12)')  &
                        'in cloudm: DCW1,DQRW3= ',dcw1, dqrw3
      endif

!---------------------------------------------------------------------


end subroutine don_cm_move_parcel_k




!######################################################################

subroutine don_cm_lcl_k    &
         (Param, t_init, p_init, mr_init, t_lcl, p_lcl, mr_lcl,  &
          lcl_reached, ermesg, error)

!---------------------------------------------------------------------
!    subroutine don_cm_lcl_k computes the lifting conden-
!    sation level by raising a parcel with temperature t_init and vapor 
!    mixing ratio mr_init adiabatically from pressure p_init until satur-
!    ation is reached at pressure p_lcl, temperature t_lcl and vapor 
!    mixing ratio mr_lcl. a flag lcl_reached is set to .true. to indicate
!    that the lcl was successfully reached.
!---------------------------------------------------------------------

use donner_types_mod, only : donner_param_type
use sat_vapor_pres_k_mod, only: compute_mrs_k

implicit none

!-----------------------------------------------------------------------
type(donner_param_type), intent(in)   :: Param
real,                    intent(in)   :: t_init, p_init, mr_init
real,                    intent(out)  :: t_lcl, p_lcl, mr_lcl
logical,                 intent(out)  :: lcl_reached
character(len=*),        intent(out)  :: ermesg
integer,                 intent(out)  :: error

!---------------------------------------------------------------------
!   intent(in) variables:
!
!      t_init         temperature of parcel at its launch point
!                     [ deg K ]
!      p_init         pressure of the parcel at its launch point [ Pa ]
!      mr_init        mixing ratio of the parcel at its launch point
!                     [ kg(h2o) / kg(air) ]
!
!   intent(out) variables:
!
!      t_lcl          temperature at lifting condensation level 
!                     [ deg K ]
!      p_lcl          pressure at lifting condensation level [ Pa ]
!      mr_lcl         mixing ratio of parcel at lifting condensation
!                     level [ kg(h2o) / kg(air) ]
!      lcl_reached    lcl was reached for this parcel ?
!
!----------------------------------------------------------------------


!---------------------------------------------------------------------
!   local variables:

      real    ::  t_parcel, p_start, p_end, gamma, rs, kappa_moist
      integer ::  max_levels
      integer ::  k, nbad

!---------------------------------------------------------------------
!   local variables:
!
!      t_parcel     temperature of parcel at current level
!      p_start      pressure at start of current iteration [ Pa ]
!      p_end        pressure at end of current iteration [ Pa ]
!      gamma        adiabatic lapse rate for moist air at current 
!                   location [ deg K / Pa ] 
!      es           saturation vapor pressure at temperature t_parcel 
!                   [ Pa ]
!      rs           saturation mixing ratio at t_parcel and p_end 
!                   [ kg(h2o) / kg(air) ]
!      kappa_moist  exponent in expression for moist potential 
!                   temperature [ dimensionless ]
!      max_levels   maximum number of iterations of parcel movement
!                   which must be taken before knowing that the
!                   parcel will not undergo deep convection
!      k            do-loop index
!
!----------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' '  ; error = 0

!---------------------------------------------------------------------
!    define an approximation for kappa (R/Cp) for moist air 
!    (kappa_moist). 
!---------------------------------------------------------------------
      kappa_moist = (Param%RDGAS/Param%CP_AIR)*(1. +    &
                    (Param%RVGAS/Param%RDGAS -     &
                     Param%CP_VAPOR/Param%CP_AIR)*mr_init)
      
!---------------------------------------------------------------------
!    define the initial values for t_parcel (the parcel temperature) 
!    and p_start and p_end, which define the pressure increment through
!    which the parcel rises during the first iteration.
!---------------------------------------------------------------------
      t_parcel = t_init
      p_start  = p_init
      p_end    = p_start + Param%parcel_dp

!---------------------------------------------------------------------
!    initialize the output variables.
!---------------------------------------------------------------------
      p_lcl  = 0.
      t_lcl  = 0.
      mr_lcl = 0.
      lcl_reached = .false.

!---------------------------------------------------------------------
!    define the maximum number of increments (max_levels) that the 
!    parcel may be raised without reaching condensation before it is 
!    no longer viable as a deep convection parcel.
!---------------------------------------------------------------------
      max_levels = int( (p_init - Param%upper_limit_for_lfc)/  &
                        abs(Param%parcel_dp) ) + 1

!----------------------------------------------------------------------
!    raise the parcel in increments of parcel_dp, following the adia-
!    batic lapse rate for moist air, until either the lcl is reached
!    or the parcel is no longer capable of undergoing deep convection. 
!----------------------------------------------------------------------
      do k=1,max_levels

!---------------------------------------------------------------------
!    exit the loop if the parcel has risen beyond where deep convection
!    is possible.
!---------------------------------------------------------------------
        if (p_end < Param%upper_limit_for_lfc) exit

!---------------------------------------------------------------------
!    define the lapse rate of temp with respect to pressure for moist
!    air (gamma). update the parcel temperature to the value at the
!    top of this layer.
!---------------------------------------------------------------------
        gamma = kappa_moist*t_parcel/p_start
        t_parcel = t_parcel + gamma*Param%parcel_dp 

!---------------------------------------------------------------------
!    determine the saturation mixing ratio for this temperature and
!    pressure. 
!---------------------------------------------------------------------
        call compute_mrs_k (t_parcel, p_end,                  &
                            Param%d622 , Param%d608 , rs, nbad, &
                            mr = mr_init)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (nbad /= 0) then
          ermesg = 'subroutine don_cm_lcl_k: '// &
                   'temperatures out of range of esat table'
          error = 1
          return
        endif

!---------------------------------------------------------------------
!   determine if the parcel is saturated. if it is, define the lcl 
!   values of temperature (t_lcl), mixing ratio (mr_lcl) and pressure 
!   (p_lcl), and return to the calling routine.
!---------------------------------------------------------------------
        if (rs <= mr_init) then
          p_lcl = p_end
          t_lcl = t_parcel
          mr_lcl = mr_init                     
          lcl_reached = .true.
          exit
        endif

!---------------------------------------------------------------------
!    define the starting and ending pressures for the next iteration.
!---------------------------------------------------------------------
        p_start = p_end
        p_end   = p_end + Param%parcel_dp
      end do

!---------------------------------------------------------------------


 end subroutine don_cm_lcl_k  


!#####################################################################


subroutine don_cm_mesub_k     &
         (Nml, pfull_c, nlev_lsm, me, diag_unit, debug_ijt, Param, cu,           &
          ci_liq_cond, ci_ice_cond, pmelt_lsm, cell_precip, &
          dint, plzb_c, pb, pt_kou, temp_c, phalf_c,   &
          ca_liq, ca_ice, ecd, ecd_liq, ecd_ice, ecei_liq,   &
          ece, ece_liq, ece_ice, meso_freeze, meso_melt, ermesg, error)

!----------------------------------------------------------------------
!    subroutine mesub calculates mesoscale heat and moisture sources,
!    using a variation on the Leary and Houze (JAS, 1980) procedure.
!    the defined fields are condensate transferred from cell to anvil
!    (ca), condensate evaporated in convective downdrafts (ecd), conden-
!    sate evaporated in convective updrafts (ece), the condensate 
!    entering the anvil which has not yet been frozen (meso_freeze), and 
!    the amount of condensate which must be melted in the mesoscale down-
!    draft to assure ice conservation (meso_melt). the subroutine is 
!    called separately for each ensemble member. for notation, see 
!    "Cu Closure A notes," 2/97.
!----------------------------------------------------------------------

use donner_types_mod, only : donner_param_type, donner_nml_type

implicit none

!----------------------------------------------------------------------
type(donner_nml_type),         intent(in)    :: Nml
integer,                       intent(in)    :: nlev_lsm, me, diag_unit
logical,                       intent(in)    :: debug_ijt
type(donner_param_type),       intent(in)    :: Param
real,                          intent(in)    :: cu, cell_precip, dint, &
                                                plzb_c, pb, pt_kou
real,   dimension(nlev_lsm),   intent(in)    :: temp_c, pfull_c
real,   dimension(nlev_lsm+1), intent(in)    :: phalf_c
real,                          intent(out)   :: ca_liq, ca_ice
real,                          intent(in)    :: pmelt_lsm, &
                                                ci_liq_cond, ci_ice_cond
real,   dimension(nlev_lsm),   intent(out)   :: ecd, ece, meso_freeze, &
                                                meso_melt, &
                                                ecd_liq, ecd_ice, &
                                                ece_liq, ece_ice
real,                          intent(out)   :: ecei_liq
character(len=*),              intent(out)   :: ermesg
integer,                       intent(out)   :: error

!---------------------------------------------------------------------
!   intent(in) variables:
!
!       cu           column integrated condensation integral
!                    [ mm / day ]
!       cell_precip  column integrated precipitation integral
!                    [ mm / day ]
!       dint???      water mass frozen in convective updraft
!            ??????  plus ice deposited convective updraft
!                    [ kg(h2o) /( (m**2) sec) ]
!                    weighted as cu,cell_precip
!       plzb_c       pressure at level of zero buoyancy [ Pa ]
!       ps           surface pressure [ Pa ]
!       pb           cloud-base pressure [ Pa ]
!       pt_kou       cloud-top pressure [ Pa ]
!       pmelt_lsm    pressure at bottom of layer in which melting 
!                    begins   [ Pa ]
!       phalf_c      large-scale model pressure half-levels (Pa)
!       debug_ijt    is this a diagnostics column ?
!       diag_unit    output unit number for this diagnostics column
!
!   intent(out) variables:
!
!       ca           total condensate transfered from cells to anvil 
!                    by this ensemble member [ mm/day ]
!       ecd          profile of condensate evaporated in convective
!                    downdraft on large-scale model grid 
!                    [ g(h2o) / kg(air) / day ] 
!       ece          profile of condensate evaporated in convective 
!                    updraft on large-scale model grid 
!                    [ g(h2o) / kg(air) / day ] 
!       meso_freeze  profile of condensate which is frozen upon enter-
!                    ing the anvil on the large-scale grid
!                    [ g(h2o) / kg(air) / day ] 
!       meso_melt    profile of condensate which is melted in mesoscale
!                    downdraft on large-scale model grid
!                    [ g(h2o) / kg(air) / day ] 
!
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!   local variables:
 
      integer ::     k
      real    ::  avail_meso_cd     ! fraction of column integrated
                                    ! condensation available to meso-
                                    ! scale circulation (1. - gnu)
                                    ! [ dimensionless ]
      real    ::  caa               ! amount of condensate which must
                                    ! be frozen when it enters the anvil
                                    ! [ g(h2o) / kg(air) / day ]
      real    ::  dint2             ! amount of condensate which has
                                    ! been frozen in the cumulus updraft
                                    ! before entering the anvil
                                    ! [ g(h2o) / kg(air) / day ]
      real    ::  ecda              ! amount of condensate evaporated 
                                    ! in cumulus downdrafts
                                    ! [ g(h2o) / kg(air) / day ]
      real :: ecda_liq, ecda_ice
      real    ::  ecdi              ! amount of condensate evaporated 
                                    ! in cumulus downdrafts [ mm / day ]
      real :: ecdi_liq, ecdi_ice
      real    ::  ecea              ! amount of condensate evaporated 
                                    ! in cumulus updrafts 
                                    ! [ g(h2o) / kg(air) / day ]
      real :: ecea_liq, ecea_ice
      real    ::  ecei              ! amount of condensate evaporated 
                                    ! in cumulus updrafts [ mm / day ]
      real ::           ecei_ice
      real    ::  elta              ! amount of condensate which must
                                    ! be melted in the mesoscale down-
                                    ! draft to conserve ice mass
                                    ! [ g(h2o) / kg(air) / day ]
      real    ::  gnu               ! fraction of column integrated 
                                    ! condensation which precipitates
                                    ! out [ dimensionless ]
      real    ::  ptt               ! pressure one cloud model delta p 
                                    ! above cloud top [ Pa ]
      real    ::  pzm               ! pressure at base of mesoscale 
                                    ! circulation [ Pa ]
      real    ::  pztm              ! pressure at top of mesoscale cir-
                                    ! culation [ Pa ]
      real    ::  p1                ! lower pressure limit for the layer
                                    ! in which one of the physical
                                    ! processes is occurring [ Pa ]
      real    ::  p2                ! upper pressure limit for the layer
                                    ! in which one of the physical
                                    ! processes is occurring [ Pa ]
      integer :: itrop
      real    :: ptrop

!---------------------------------------------------------------------
!   local variables:
!
!      
      ermesg = '  ' ; error = 0

!---------------------------------------------------------------------
!    define pressure one cloud-model level above cloud top (ptt). 
!    define the pressure at top of mesoscale updraft (pztm, 300 hPa 
!    plus one model-layer pressure thickness above cloud top).
!---------------------------------------------------------------------
      ptt = pt_kou + Param%dp_of_cloud_model
      pztm = ptt - 300.E02

!---------------------------------------------------------------------
!    restrict pztm to >= 100 hPa, cf Ackerman et al (JAS,1988), unless 
!    pt_kou <= 100 hPa. it was found in AM2p9 that the stratospheric 
!    water vapor was excessive with this pztm restriction, so pztm is now
!    set to be no higher than the level of zero buoyancy, or if the
!    cloud top is above the level of zero buoyancy, it is set to one 
!    model layer above the level of zero buoyancy. 
!---------------------------------------------------------------------
      if (pztm < plzb_c) pztm = plzb_c
      if (ptt < plzb_c)  pztm = plzb_c + Param%dp_of_cloud_model

      if (Nml%limit_pztm_to_tropo) then
        call find_tropopause (nlev_lsm, temp_c, pfull_c, ptrop, itrop)
        pztm = MAX (pztm, ptrop)
      endif

!---------------------------------------------------------------------
!    define the base of the mesoscale updraft (pzm), as the layer imm-
!    ediately above cloud top, or, if the top of the mesoscale updraft
!    has been redefined to be at or just above the level of zero 
!    buoyancy, to be one layer below the mesoscale updraft top. 
!---------------------------------------------------------------------
      pzm = ptt
      if (pzm <= pztm) pzm = pztm - Param%dp_of_cloud_model

!---------------------------------------------------------------------
!    if in a diagnostics column, output the convective rain 
!    (cell_precip), convective updraft condensation (cu), and the pres-
!    sure at the level of zero buoyancy (plzb_c).
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)') 'in mesub: rc,cu= ',  &
                                                    cell_precip, cu
        write (diag_unit, '(a,  e20.12)') 'in mesub: plzb = ',plzb_c
      endif

!----------------------------------------------------------------------
!    define the ratio of precipitation to condensation for the current
!    ensemble member (gnu). define the remaining fraction of condens-
!    ation 1 - gnu as the condensate available to the mesoscale circ-
!    ulation (avail_meso_cd). define the mass of this available conden-
!    sate which is evaporated in convective downdrafts (ecdi), the mass
!    evaporated into the cell environment (ecei) and the portion incor-
!    porated into the mesoscale region (ca). this partitioning is 
!    defined by the parameters evap_in_downdraft, evap_in_environ and 
!    entrained_into_meso, taken from the work of Leary and Houze 
!    (JAS, 1980).
!----------------------------------------------------------------------
      gnu = cell_precip/cu
      avail_meso_cd = 1. - gnu
      ecdi  = (Param%evap_in_downdrafts*avail_meso_cd)*cu
      ecdi_liq  = (Param%evap_in_downdrafts*avail_meso_cd)* &
                         (Param%seconds_per_day*ci_liq_cond)
      ecdi_ice  = (Param%evap_in_downdrafts*avail_meso_cd)*     &
                         (Param%seconds_per_day*ci_ice_cond)
      ecei  = (Param%evap_in_environ*avail_meso_cd)*cu
      ecei_liq  = (Param%evap_in_environ*avail_meso_cd)*    &
                         (Param%seconds_per_day*ci_liq_cond)
      ecei_ice  = (Param%evap_in_environ*avail_meso_cd)*     &
                         (Param%seconds_per_day*ci_ice_cond)
      ca_liq    = (Param%entrained_into_meso*avail_meso_cd)*   &
                         (Param%seconds_per_day*ci_liq_cond)
      ca_ice    = (Param%entrained_into_meso*avail_meso_cd)*   &
                         (Param%seconds_per_day*ci_ice_cond)
      if (debug_ijt) then
        write (diag_unit, '(a, 3e20.12)')  &
               'in mesub: cu, h1_liqintg, h1_iceintg= ',  &
            cu, ci_liq_cond*Param%seconds_per_day      ,  &
                 ci_ice_cond*Param%seconds_per_day      
      endif

!---------------------------------------------------------------------
!    if in a diagnostics column, output the ratio of convective rain 
!    to convective updraft condensation (gnu) and the mass entrained
!    into the mesoscale region (ca).
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12)')  'in mesub: gnu= ',gnu
        write (diag_unit, '(a, e20.12)') 'in mesub: ca= ',  &
                                                       ca_liq + ca_ice
        write (diag_unit, '(a, 2e20.12)') 'in mesub: ca_liq,ca_ice= ', &
                                                ca_liq, ca_ice
      endif

!--------------------------------------------------------------------
!    calculate the mass of water which must be frozen as it enters the
!    mesoscale anvil (caa). if no freezing has occurred in the cumulus
!    updraft (i.e., dint2 = 0) then this will be ca, the total mass 
!    available to the anvil. if freezing has occurred, (ie, 
!    dint2 /= 0.), then the amount to be frozen is the total amount 
!    available (ca) plus additional vapor mass deposited on the ice in 
!    the updraft (ecei), less that which has already frozen (dints). 
!    dints and caa are expressed in units of g(h2o) per kg(air) per day.
!--------------------------------------------------------------------
!9/15/07, 1037AM:
      dint2 = avail_meso_cd*(dint                               )*  &
               8.64e07*Param%grav/(pzm - pztm)

      if (dint2 /= 0.)  then 
       caa = ((ca_liq + ecei_liq)*Param%grav*1000./(pzm - pztm)) - dint2
      else
        caa = ca_liq*Param%grav*1000./(pzm - pztm)
      endif

!---------------------------------------------------------------------
!    if in a diagnostics column, output the previously frozen condensate
!    (dint2), the additional amount to be frozen (caa) and the pressure
!    range over which the freezing will occur (pzm, pztm). if 
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a,  e20.12)')  &
                         'in mesub:     dint           =',    dint 
        write (diag_unit, '(a, 2e20.12)')  &
                         'in mesub:     dint2, ecei_liq=',    dint2, &
                                                          ecei_liq
        write (diag_unit, '(a, 3e20.12)')  &
                           'in mesub: caa,pzm,pztm= ',caa,pzm,pztm
      endif

!---------------------------------------------------------------------
!    if there is additional condensate which must be frozen upon enter-
!    ing the anvil, call map_hi_res_intgl_to_lo_res_col to spread this 
!    additional freezing uniformly over the region between anvil base 
!    (pzm) and anvil top (pztm) in the large-scale model. store the out-
!    put in array meso_freeze. if no additional freezing is needed, set 
!    meso_freeze to be 0.0.
!---------------------------------------------------------------------
      if (caa > 0.)  then 
        if (debug_ijt) then
          write (diag_unit, '(a, e20.12, 2f19.10)')  &
                      'in cm_intgl_to_gcm_col: xav,p1,p2= ',caa, pzm, &
                                                  pztm
        endif
        call don_u_map_hires_i_to_lores_c_k   &
             (nlev_lsm, caa, pzm, pztm, phalf_c, meso_freeze, ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (error /= 0 ) return
        if (debug_ijt) then
          do k=1,nlev_lsm       
            if (meso_freeze(k) /= 0.0) then
              write (diag_unit, '(a, i4, e20.12)') &
                    'in cm_intgl_to_gcm_col: k,x= ',k, meso_freeze   (k)
            endif
          end do
        endif
      else
        meso_freeze = 0.
      endif

!---------------------------------------------------------------------
!    define the evaporation which occurs in the convective downdraft.
!    the convective downdraft is assumed to originate one layer above
!    the cloud top (ptt) and extend to the surface (phalf_c(1)). 
!    convert the convective downdraft evaporation to units of
!    g(h20) / kg(air) per day.
!---------------------------------------------------------------------
      ecda = ecdi*Param%grav*1000./(phalf_c(1) - ptt)
      ecda_liq = ecdi_liq*Param%grav*1000./(phalf_c(1) - ptt)
      ecda_ice = ecdi_ice*Param%grav*1000./(phalf_c(1) - ptt)

!---------------------------------------------------------------------
!    if in a diagnostics column, output the convective downdraft evap-
!    oration (ecda) and the large-scale model pressure limits over which
!    this evaporation occurs (phalf_c(1), ptt).
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 3e20.12)')  &
                         'in mesub: ecda,p1,pz0= ',ecda,phalf_c(1),ptt
        write (diag_unit, '(a, 2e20.12)')  &
                         'in mesub: ecda_liq, ecda_ice= ',  &
                                ecda_liq, ecda_ice
      endif

!---------------------------------------------------------------------
!    call map_hi_res_intgl_to_lo_res_col to spread the integrated evap-
!    oration in convective downdrafts uniformly over the region between
!    the surface (phalf_c(1)) and the anvil base (pzm) and the top of 
!    cloud (ptt). output field is ecd.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12, 2f19.10)')  &
             'in cm_intgl_to_gcm_col: xav,p1,p2= ',ecda, phalf_c(1) , &
                                                  ptt 
      endif
      call don_u_map_hires_i_to_lores_c_k   &
           (nlev_lsm, ecda, phalf_c(1), ptt, phalf_c, ecd, ermesg, error)
      call don_u_map_hires_i_to_lores_c_k   &
         (nlev_lsm, ecda_liq, phalf_c(1), ptt, phalf_c, ecd_liq, ermesg, error)
      call don_u_map_hires_i_to_lores_c_k   &
         (nlev_lsm, ecda_ice, phalf_c(1), ptt, phalf_c, ecd_ice, ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return
      if (debug_ijt) then
        do k=1,nlev_lsm       
          if (ecd(k) /= 0.0) then
            write (diag_unit, '(a, i4, e20.12)') &
                'in cm_intgl_to_gcm_col: k,ecd= ',k, ecd   (k)
            write (diag_unit, '(a, i4, 2e20.12)') &
              'in cm_intgl_to_gcm_col: k,ecdliq,ecdice= ',k, &
                       ecd_liq(k), ecd_ice(k)
          endif
        end do
      endif

!---------------------------------------------------------------------
!    be sure that the melting level in the large-scale model (pmelt_lsm)
!    is below the top of the mesoscale circulation (pztm),and above
!    cloud base (pb). if not, no melting will occur; set p2 to be 0.0.
!---------------------------------------------------------------------
      elta = 0.
      if (pmelt_lsm  < pztm                    )  then
        meso_melt = 0.
      if (debug_ijt) then
        write (diag_unit, '(a, 2f19.10)') &
                 ' NO MELTING DONE: melting level above top of &
                    &mesoscale circulation : pmelt_lsm,pztm',     &
                                                pmelt_lsm, pztm      
      endif

!---------------------------------------------------------------------
!    if pmelt_lsm is within the region of the cloud and mesoscale circ-
!    ulation, calculate any melting that must occur in the mesoscale
!    downdraft in order to conserve ice mass; ie, if the amount to be
!    frozen was calculated as more than the available condensate, then
!    the excess must be melted, and is done so in the mesoscale down-
!    draft between the melting level and cloud base.
!---------------------------------------------------------------------
      else if (pmelt_lsm >= pztm .and. pmelt_lsm <= pb) then
        p2 = pmelt_lsm
        p1 = pb
        if (caa <= 0.) then 
          caa = -caa*(pzm - pztm)/(pb - p2)
          elta = caa
        endif
      if (debug_ijt) then
        write (diag_unit, '(a, 3f19.10)') &
                   'MELTING DONE: pmelt_lsm,pb,caa  ',pmelt_lsm, pb, &
                                                            caa  
      endif

!---------------------------------------------------------------------
!    if in diagnostics column, output the melting (elta) and the 
!    pressures defining the layer in which it occurs (pb, p2)
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12, 2f19.10)') &
                           'in mesub: elta,p1,p2= ',elta,p1,p2
      endif

!---------------------------------------------------------------------
!    call map_hi_res_intgl_to_lo_res_col to spread the required melting
!    resulting from excessive freezing over the layer between cloud base
!    and the melting level. output field is meso_melt.
!---------------------------------------------------------------------
      call don_u_map_hires_i_to_lores_c_k   &
           (nlev_lsm, elta, p1, p2, phalf_c, meso_melt, ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return
      if (debug_ijt) then
        do k=1,nlev_lsm       
          if (meso_melt(k) /= 0.0) then
            write (diag_unit, '(a, i4, e20.12)') &
                 'in cm_intgl_to_gcm_col: k,meso_melt= ',k, meso_melt(k)
          endif
        end do
      endif

      else if (pmelt_lsm > pb) then
        meso_melt = 0.
        if (pmelt_lsm == phalf_c(1)) then
          if (debug_ijt) then
             write (diag_unit, '(a)') &
                   'NO MELTING LEVEL PRESENT IN COLUMN'
          endif
        else
! melt below cloud base 
      if (debug_ijt) then
        write (diag_unit, '(a, 2f19.10)') &
            ' NO MELTING DONE: melting level below PB: pmelt_lsm,pb', &
                                                      pmelt_lsm, pb
      endif
      endif
      endif ! (pmelt<pztm or pmelt > pb)

!---------------------------------------------------------------------
!    calculate the evaporation which occurs in the convective 
!    updraft.
!    this is spread between 50 hPa below cloud top and 10 hPa above 
!    cloud top.
!---------------------------------------------------------------------
      p1 = pt_kou + 50.0e02
      p2 = ptt
      ecea = ecei*Param%grav*1000./(p1-p2)
      ecea_liq = ecei_liq*Param%grav*1000./(p1-p2)
      ecea_ice = ecei_ice*Param%grav*1000./(p1-p2)

!---------------------------------------------------------------------
!    if in diagnostics column, output the convective updraft evaporation
!    (ecea, ecei) and the large-scale model pressure layer limits over 
!    which it occurs (p1, p2).
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)')  &
                         'in mesub: ecea,ecei= ',ecea, ecei
        write (diag_unit, '(a, 2e20.12)')  &
                         'in mesub: LIQecea,ecei= ',ecea_liq, ecei_liq
        write (diag_unit, '(a, 2e20.12)')  &
                         'in mesub: ICEecea,ecei= ',ecea_ice, ecei_ice
        write (diag_unit, '(a, e20.12, 2f19.10)')  &
                         'in mesub: ecea,p1,p2= ',ecea, p1, p2
      endif

!---------------------------------------------------------------------
!    call map_hi_res_intgl_to_lo_res_col to spread the integrated evap-
!    oration in convective updrafts uniformly over the designated 
!    region.  output field is ece.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12, 2f19.10)')  &
                      'in cm_intgl_to_gcm_col: xav,p1,p2= ',ecea, p1, &
                                                  p2
      endif
      call don_u_map_hires_i_to_lores_c_k   &
           (nlev_lsm, ecea, p1, p2, phalf_c, ece, ermesg, error)
      call don_u_map_hires_i_to_lores_c_k   &
           (nlev_lsm, ecea_liq, p1, p2, phalf_c, ece_liq, ermesg, error)
      call don_u_map_hires_i_to_lores_c_k   &
           (nlev_lsm, ecea_ice, p1, p2, phalf_c, ece_ice, ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return
      if (debug_ijt) then
        do k=1,nlev_lsm     
          if (ece(k) /= 0.0) then
            write (diag_unit, '(a, i4, e20.12)') &
                           'in cm_intgl_to_gcm_col: k,x= ',k, ece   (k)
          endif
        end do
      endif

!---------------------------------------------------------------------


end subroutine don_cm_mesub_k



!######################################################################

subroutine don_cm_compute_vert_fluxes_k   &
         (nlev_hires, ntr, ncc_kou, kou, diag_unit, debug_ijt, &
          do_donner_tracer, Param, pt_kou, cld_press, rcl, te, mre, wv, &
          tcc, dpf, xclo, xtrae, apt, efchr, emfhr, etfhr, ermesg, error)

!---------------------------------------------------------------------
!    subroutine compute_vertical_fluxes computes vertical flux conver-
!    gence terms and their column integrals for temperature, entropy, 
!    moisture and tracers within the cloud. 
!---------------------------------------------------------------------

use donner_types_mod, only : donner_param_type
use sat_vapor_pres_k_mod, only: compute_qs_k

implicit none 

!---------------------------------------------------------------------
integer,                            intent(in)     :: nlev_hires, ntr,  &
                                                      ncc_kou, kou,  &
                                                      diag_unit
logical,                            intent(in)     :: debug_ijt,   &
                                                      do_donner_tracer
type(donner_param_type),            intent(in)     :: Param
real,                               intent(in)     :: pt_kou
real,    dimension(nlev_hires),     intent(in)     :: cld_press, rcl,  &
                                                      te, mre, wv, tcc, &
                                                      dpf
real,    dimension(nlev_hires,ntr), intent(in)     :: xclo, xtrae      
real,                               intent(out)    :: apt               
real,    dimension(nlev_hires),     intent(out)    :: efchr, emfhr    
real,    dimension(nlev_hires,ntr), intent(out)    :: etfhr            
character(len=*),                   intent(out)    :: ermesg
integer,                            intent(out)    :: error

!---------------------------------------------------------------------
!   intent(in) variables:
!
!      pt_kou      pressure at cloud top [ Pa ]
!      cld_press   pressure at full levels in cloud model [ Pa ]
!      rcl         vertical profile of cloud radius [ m ]
!      te          environmental temperature profile [ deg K ]
!      mre         environmental vapor mixing ratio profile 
!                  [ kg(h2o) / kg (dry air) ]
!      wv          in-cloud vertical velocity profile [ m / sec ]
!      tcc         in-cloud temperature profile [ deg K ] 
!      xclo        in-cloud tracer profiles 
!                  [ kg(tracer) / kg (dry air) ]
!      xtrae       environmental tracer profiles 
!                  [ kg(tracer) / kg (dry air) ]
!      debug_ijt   logical indicating whether diagnostics are desired
!                  for column
!      ncc_kou     cloud model level index of level at or above 
!                  cloud top
!      kou         index of current ensemble member 
!      diag_unit   output unit for diagnostics for this column
!
!   intent(out) variables:
!
!      apt         ratio of cloud area at cloud top to that at cloud 
!                  base [ dimensionless ]
!      efchr       vertical entropy flux convergence [ deg K / sec ]
!      emfhr       vertical moisture flux convergence 
!                  [ kg(h2o) / ( kg(dry air) sec ] 
!      etfhr       vertical tracer flux convergence 
!                  [ kg(tracer) / ( kg(dry air) sec ] 
!
!----------------------------------------------------------------------

!---------------------------------------------------------------------
!   local variables:

      real, dimension(nlev_hires)       ::  dpdz, q_sat, ehf, emf
      real, dimension(nlev_hires,ntr) ::  etf                      
      real, dimension(ntr)      ::  sumetf
      real                        ::  p, pl, ph, dpp, ptt, exf, &
                                      thetf, tv_env, tv_cld,  &
                                      dpdz_cb, sumemf, sumefc, sumthet
      integer                     ::  kcl, kch
      integer                     ::  kc, kcont, nbad


!---------------------------------------------------------------------
!   local variables:
!
!        dpdz            value of dp/dz, evaluated from hydrostatic
!                        equation, with virtual mass coefficient to
!                        account for non-hydrostatic effects. see
!                        Donner, 1986, j. atmos. sci., 43, p 2288.
!                        [ kg(air) / (sec**2 m**2) ]
!        q_sat           saturation specific humidity 
!                        [ kg(h2o) / kg (air) ]
!        ehf             vertical temperature flux 
!                        [ ( kg(air) deg K ) / (m sec**3) ]
!        emf             vertical moisture flux
!                        [ kg (h2o) / (m sec**3) ]
!        etf             vertical tracer flux
!                        [  kg(tracer) / (m sec**3) ]
!        sumetf          column integral of vertical tracer flux 
!                        convergence [  kg(tracer) / (m sec**3) ]
!        p               pressure  at level where subgrid scale temper-
!                        ature flux is calculated [ Pa ]
!        pl              pressure at lower flux interface [ Pa ]
!        ph              pressure at upper flux interface [ Pa ]
!        dpp             pressure depth over which flux is being 
!                        calculated [ Pa ]
!        ptt             pressure one level above cloud top [ Pa ]
!        exf             perturbation vertical subgrid scale temp-
!                        erature flux [ ( kg(air) deg K ) / (m sec**3) ]
!        thetf           vertical flux of potential temperature at
!                        cloud base [ ( kg(air) deg K ) / (m sec**3) ]
!        esat            saturation vapor pressure [ Pa ]
!        tv_env          environmental virtual temperature [ deg K ]
!        tv_clda         in-cloud virtual temperature [ deg K ]
!        dpdz_cb         value of dp/dz over the lowest cloud model
!                        pressure layer [ kg(air) / (sec**2 m**2) ]
!        sumemf          column integral of vertical moisture flux 
!                        convergence [ kg (h2o) / (m sec**3) ]
!        sumefc          column integral of vertical entropy flux 
!                        convergence [ ( kg(air) deg K ) / (m sec**3) ]
!        sumthet         column integral of vertical potential temper-
!                        ature flux convergence 
!                        [ ( kg(air) deg K ) / (m sec**3) ]
!        kcl             k index of lower interface level
!        kch             k index of upper interface level
!        kc, kcont       do-loop indices
!
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!--------------------------------------------------------------------
!    initialize the output arrays.
!--------------------------------------------------------------------
      do kc=1,nlev_hires
        do kcont=1,ntr
          etfhr(kc,kcont) = 0.
        end do
        efchr(kc) = 0.
        emfhr(kc) = 0.
      end do

!---------------------------------------------------------------------
!    initialize variables which will collect column sums.
!---------------------------------------------------------------------
      sumemf    = 0.
      sumefc    = 0.
      sumthet   = 0.
      sumetf(1:ntr) = 0.

!--------------------------------------------------------------------
!    initialize the ratio of cloud area at cloud top to cloud area at 
!    cloud base. 
!--------------------------------------------------------------------
      apt = 0.

!---------------------------------------------------------------------
!    loop over the cloud model levels from cloud base (kc=1) to the 
!    level just below cloud top (ncc_kou-1), defining variables needed 
!    for the flux convergence calculations.
!---------------------------------------------------------------------
      do kc=1,ncc_kou-1 

!---------------------------------------------------------------------
!    define the environmental and in-cloud virtual temperatures. NOTE: 
!    q_sat is correctly a specific humidity; mre(kc) is incorrectly a 
!    mixing ratio.
!---------------------------------------------------------------------
        call compute_qs_k (tcc(kc), cld_press(kc), Param%D622,  &
                           Param%D608, q_sat(kc), nbad)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (nbad /= 0) then
          ermesg = 'subroutine don_cm_compute_vert_fluxes_k: '// &
                   'temperatures out of range of esat table'
          error = 1
          return
        endif
        tv_env     = te(kc)* (1. + Param%d608*(mre(kc)/(1. + mre(kc))))
        tv_cld     = tcc(kc)*(1. + Param%d608*q_sat(kc))

!---------------------------------------------------------------------
!    compute dp/dz (dpdz) as in equation (B3) in Donner (1986), 
!    j. atm. sci., 43, pp 2277-2288.  here the expression uses the 
!    pressure mid-way between cloud base and level 2. compute vertical 
!    eddy transport of temperature (ehf), moisture (emf) and tracers
!    (etf), using equation (3) from Donner et al (1982), j. atm. sci.,
!    39, pp 2159-2181. multiply by dpdz to put flux in pressure (omega) 
!    units. note (1 - a) is approximated as 1. note that q_sat is a 
!    specific humidity while mre is a mixing ratio.
!---------------------------------------------------------------------
        if (kc == 1) then

!--------------------------------------------------------------------
!    at cloud base level, calculate dpdz valid over the interval between
!    level 1 and level 2 (dpdz_cb). this term is needed to determine 
!    the subgrid scale vertical temperature flux valid in the first 
!    cloud layer.
!--------------------------------------------------------------------
          dpdz_cb = -Param%grav*(0.5*(cld_press(1) + cld_press(2)))* &
                    (Param%virt_mass_co*tv_env + tv_cld)/(Param%rdgas*  &
                    (1. + Param%virt_mass_co)*tv_env*tv_cld)

!----------------------------------------------------------------------
!    when at cloud top, calculate dpdz and the cloud perturbation flux 
!    terms that are valid for the model layer including cloud top. save
!    these values as index ncc_kou of the various arrays.
!--------------------------------------------------------------------
        else if (kc == ncc_kou-1) then
          dpdz(ncc_kou) = -Param%grav*pt_kou*(    &
                          Param%virt_mass_co*tv_env + tv_cld)/   &
                          (Param%rdgas*(1. + Param%virt_mass_co)*  &
                          tv_env*tv_cld)
          ehf(ncc_kou) = ((rcl(ncc_kou-1)/Param%cloud_base_radius)**2)* &
                         wv(ncc_kou-1)*dpdz(ncc_kou)*   &
                         (tcc(ncc_kou-1) - te(ncc_kou-1))
          emf(ncc_kou) = ((rcl(ncc_kou-1)/Param%cloud_base_radius)**2)* &
                         wv(ncc_kou-1)*dpdz(ncc_kou)*    &
                         (q_sat(ncc_kou-1)/(1. - q_sat(ncc_kou-1)) - &
                                                         mre(ncc_kou-1))
          do kcont=1,ntr
            etf(ncc_kou,kcont) = ((rcl(ncc_kou-1)/  &
                                 Param%cloud_base_radius)**2)* &
                                 wv(ncc_kou-1)*dpdz(ncc_kou)*  &
                                 (xclo(ncc_kou-1,kcont) -   &
                                                  xtrae(ncc_kou-1,kcont))
          end do
        endif

!---------------------------------------------------------------------
!    compute dpdz and the vertical eddy flux terms at interior cloud
!    levels.
!---------------------------------------------------------------------
        dpdz(kc) = -Param%grav*cld_press(kc)*  &
                   (Param%virt_mass_co*tv_env + tv_cld)/&
                   (Param%rdgas*(1. + Param%virt_mass_co)*tv_env*tv_cld)
        ehf(kc) = ((rcl(kc)/Param%cloud_base_radius)**2)*wv(kc)*  &
                  dpdz(kc)*(tcc(kc) - te(kc))
        emf(kc) = ((rcl(kc)/Param%cloud_base_radius)**2)*wv(kc)*   &
                  dpdz(kc)*(q_sat(kc)/(1. - q_sat(kc)) - mre(kc))
        do kcont=1,ntr
          etf(kc,kcont) = ((rcl(kc)/Param%cloud_base_radius)**2)*  &
                            wv(kc)*dpdz(kc)*  &
                            (xclo(kc,kcont) - xtrae(kc,kcont))
        end do
      end do

!--------------------------------------------------------------------
!    loop over levels from cloud base to the highest level within
!    the cloud to compute the eddy flux convergence profiles.
!--------------------------------------------------------------------
      do kc=1,ncc_kou-1  

!--------------------------------------------------------------------
!    define the lower interface index (kcl) to be one less than 
!    the current level. at cloud base, define the lower interface index
!    to be the current index. define the lower interface pressure (pl).
!    if the lower interface is at or above cloud top, calculation is
!    complete in this column, so return. define the upper interface 
!    pressure (ph) to be the pressure one level above the current index.
!    define the delta p between the interface pressures (dpp).
!--------------------------------------------------------------------
        kcl = MAX (1, kc-1)
        pl = cld_press(kcl)
        if (pl <= pt_kou) return
        ph = cld_press(kc+1)
        dpp = (ph - pl)/2.0

!---------------------------------------------------------------------
!     if the upper interface pressure is not above cloud top, set the 
!     upper interface index (kch) to be one above the current index. set
!     the current level pressure p to be the average value of the higher
!     and lower level pressures, which will be the current level value,
!     except at the lowest level. 
!---------------------------------------------------------------------
        if (ph >= pt_kou) then
          kch = kc + 1
          p = (pl + ph)/2.0

!---------------------------------------------------------------------
!     if the upper interface pressure is above cloud top, set the 
!     upper interface index (kch) to be the current value (kc), and set
!     the upper interface and current level pressures to be the cloud 
!     top pressure.  define the ratio of cloud-top cloud area to cloud-
!     base cloud area (apt).
!---------------------------------------------------------------------
        else
          kch = kc
          ph = pt_kou
          p = pt_kou
          apt = (rcl(kc)/rcl(1))**2
        endif

!!$$^%S^&A^^^**
!!    SHOULDN'T dpp be defined here, after any adjustnment to ph has 
!     BEEN MADE ???? Problem would be at topmost layer.
!!$$^%S^&A^^^**

!---------------------------------------------------------------------
!    if in a debug column, output values of cloud-level indices and 
!    cloud temperature, along with environmental mixing ratio and
!    temperature at these levels.
!---------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, 3i4)')  &
                             'in mulsub: kc,kcl,kch= ',kc,kcl,kch
          write (diag_unit, '(a, 3e20.14)')  &
                         'in mulsub: TCC = ',TCC(kc),TCC(kcl),TCC(kch)
          write (diag_unit, '(a, 3e20.12)')  &
                              'in mulsub: QE= ',mrE(kc),mrE(kcl),mrE(kch)
          write (diag_unit, '(a, 3e20.12)')   &
                              'in mulsub: TE= ',TE(kc),TE(kcl),TE(kch)
        endif

!----------------------------------------------------------------------
!    define the perturbation vertical subgrid scale temperature flux 
!    associated with the cloud (exf). note the special definitions at
!    cloud bottom and top.
!----------------------------------------------------------------------
        if (kc == 1) then
          exf = Param%rdgas*wv(kc)*dpdz_cb*(tcc(kc) - te(kc))*  &
                ((rcl(kc)/Param%cloud_base_radius)**2)/(Param%cp_air*p)
        else if (kc == kch) then
          exf = Param%rdgas*wv(kc)*dpdz(ncc_kou)*(tcc(kc) - te(kc))*  &
                ((rcl(kc)/Param%cloud_base_radius)**2)/(Param%cp_air*p)
        else             
          exf = Param%rdgas*wv(kc)*dpdz(kc)*(tcc(kc) - te(kc))*   &
                ((rcl(kc)/Param%cloud_base_radius)**2)/(Param%cp_air*p)
        endif

!----------------------------------------------------------------------
!    at the cloud top, define the cumulus vertical-flux convergence of 
!    temperature, moisture and tracers as the flux transport entering 
!    the layer which contains cloud top. this "through the top" trans-
!    port is assumed to be all deposited in this layer which includes 
!    cloud top. set the vertical eddy transports of temperature, moist-
!    ure and tracer at the uppermost interior cloud level to zero.
!----------------------------------------------------------------------
        if (kc == kch) then
          efchr(ncc_kou ) = ehf(ncc_kou)/Param%dp_of_cloud_model
          emfhr(ncc_kou ) = emf(ncc_kou)/Param%dp_of_cloud_model
          ehf(ncc_kou-1) = 0.
          emf(ncc_kou-1) = 0.

!##$%$%$%^&^&!!&&*** 
!!!???? is this correct ? shouldn't the flux across the upper interface
!!! of this layer be accounted for ?  I..e., 

!!!!     efchr(ncc-1) = exf +  (ehf(ncc-2) - ehf(ncc))/(2.*dpp)

!!   As it stands, it appears that all the flux into the ncc-1 layer 
!!   (ehf(ncc-2)) remains. however, some (ehf(ncc)) is exported 
!    into the ncc layer (above cloud top), suggesting a non-conservation
!!!   condition.
!##$%$%$%^&^&!!&&*** 

          efchr(ncc_kou-1) = exf + (ehf(ncc_kou-2))/(2.*dpp)
          emfhr(ncc_kou-1) = emf(ncc_kou-2)/(2.*dpp)
          do kcont=1,ntr
            etfhr(ncc_kou,kcont) = etf(ncc_kou,kcont)/Param%dp_of_cloud_model
            etf(ncc_kou-1,kcont) = 0.
            etfhr(ncc_kou-1,kcont) = etf(ncc_kou-2,kcont)/(2.*dpp)
          end do

!----------------------------------------------------------------------
!    add the appropriately pressure-weighted contributions from the
!    layer containing cloud top to the integrals of the in-cloud vert-
!    ical-flux convergence of entropy, (sumefc) potential temperature
!    (sumthet), moisture (sumemf) and tracers (sumetf).
!----------------------------------------------------------------------
          sumefc = sumefc + efchr(ncc_kou)*Param%dp_of_cloud_model/2.
          ptt = pt_kou + Param%dp_of_cloud_model
          sumthet = sumthet + efchr(ncc_kou)*    &
                                   ((1.0e05/ptt)**Param%kappa)* &
                                               Param%dp_of_cloud_model/2.
          sumemf = sumemf + emfhr(ncc_kou)*Param%dp_of_cloud_model/2.
          do kcont=1,ntr
            sumetf(kcont) = sumetf(kcont) + etfhr(ncc_kou,kcont)*  &
                            Param%dp_of_cloud_model/2.
          end do

!----------------------------------------------------------------------
!    if in diagnostic column, output the entropy and moisture flux 
!    convergence in the layer containing cloud top.
!----------------------------------------------------------------------
          if (debug_ijt) then 
            write (diag_unit, '(a, e20.12)')  &
                      'in mulsub: EFCHR(kc+1)= ',efchr(ncc_kou)
            write (diag_unit, '(a, e20.12)')  &
                          'in mulsub: EMFHRIT=kch ',emfhr(ncc_kou)
          endif

!---------------------------------------------------------------------
!    for all but the topmost layer, define the layer flux convergences
!    of entropy (efchr), moisture (emfhr) and tracers (etfhr). the ver-
!    tical flux convergence of entropy (efchr) is defined as the sum of
!    the in-cloud perturbation vertical subgrid scale temperature flux 
!    (exf) and the eddy flux convergence of temperature across 
!    the layer.
!---------------------------------------------------------------------
        else
          efchr(kc) = exf + (ehf(kcl) - ehf(kch))/(2.*dpp)
          emfhr(kc) = (emf(kcl) - emf(kch))/(2.*dpp)
          do kcont=1,ntr
            etfhr(kc,kcont) = (etf(kcl,kcont) - etf(kch,kcont))/(2.*dpp)
          end do
        endif

!----------------------------------------------------------------------
!    add the appropriately pressure-weighted contributions from the
!    current layer to the integrals of the in-cloud vertical-flux con-
!    vergence of entropy, (sumefc) potential temperature (sumthet), 
!    moisture (sumemf) and tracers (sumetf).
!----------------------------------------------------------------------
        sumefc  = sumefc  + (efchr(kc)*dpp)
        sumthet = sumthet + efchr(kc)*((1.0e05/p)**Param%kappa)*dpp
        sumemf  = sumemf  + (emfhr(kc)*dpp)
        do kcont=1,ntr
          sumetf(kcont) = sumetf(kcont) + etfhr(kc,kcont)*dpp
        end do

!---------------------------------------------------------------------
!    if in diagnostic column, write out the current vertical sum of 
!    the entropy and potential temperature flux convergence. output the
!    vertical flux of theta (thetf) and moisture (emf) at cloud base.
!    output the in-cloud (xclo) and environmental (xtrae) tracer 
!    mixing ratios.  if this is the 4th cumulus subensemble, output 
!    some fluxes, flux convergences, dpdz and other terms related to 
!    the flux calculation.
!---------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, 2e20.12)')  &
                    'in mulsub: SUMTHET,SUMEFC= ',sumthet, sumefc
          do kcont=1,ntr
            write (diag_unit, '(a, e20.12)')'in mulsub: SUMETF=', &
                                           sumetf(kcont)
          end do
          if (kc == 1) then
            thetf = ehf(kcl)*((1.0e05/pl)**Param%kappa)
            write (diag_unit, '(a, 2e20.12)')  &
                             'in mulsub: THETF,EMFF= ',thetf,emf(kcl)
          endif
          do kcont=1,ntr
            write (diag_unit, '(a, 3i4, 2e20.12)')  &
                   'in muls:kou,k,kcont,xclo,xtrae=' &
                    ,kou,kc,kcont,xclo(kc,kcont),xtrae(kc,kcont)
          end do
          if (kou == 4) then
            write (diag_unit, '(a, i4, 2f19.10)') &
                                  'in mulsub: kc,PL,PH= ',kc, pl, ph
            write (diag_unit, '(a, 4e20.12)')  &
                'in mulsub: EHFH,EHFL,EXF,efchr(kc)= ', &
                        ehf(kch), ehf(kcl), exf, efchr(kc)
            write (diag_unit, '(a, 3e20.12)')   &
                          'in mulsub: EMFH,EMFL,EMF= ',  &
                         emf(kch), emf(kcl), emfhr(kc)
            if (do_donner_tracer) then
              write (diag_unit, '(a, 3e20.12)')   &
                       'in mulsub: ETFH,ETFL,ETF=         ', &
                     etf(kch,ntr), etf(kcl,ntr), etfhr(kc,ntr)
            endif
            write (diag_unit, '(a, 3e20.12)')   &
                               'etfh diag: rcl,wv,dpdzh=         ', &
                                    rcl(kch), wv(kch), dpdz(kch)
            write (diag_unit, '(a, 3e20.12)')   &
                                'etfl diag: rcl,wv,dpdzl=         ', &
                                    rcl(kcl), wv(kcl), dpdz(kcl)
            do kcont=1,ntr
              write (diag_unit, '(a, 3e20.12)')  &
                                   'etfh diag: xclo,xtrae= ',         &
                                   xclo(kch,kcont), xtrae(kch,kcont)
              write (diag_unit, '(a, 3e20.12)')  &
                                 'etfl diag: xclo,xtrae= ',         &
                                   xclo(kcl,kcont), xtrae(kcl,kcont)
            end do
            write (diag_unit, '(a, 3e20.12)')   &
                            'in mulsub: WV,RH,QE= ',  &
                                     wv(kch), q_sat(kch), mre(kch)
            write (diag_unit, '(a, 3e20.12)')   &
                            'in mulsub: WV,RL,QE= ',   &
                                     wv(kcl), q_sat(kcl), mre(kcl)
          endif 
        endif
      end do  ! (end of kc loop)

!----------------------------------------------------------------------
!    if this is a debug column, output the moisture and temperature 
!    tendency profiles produced by this ensemble member.
!----------------------------------------------------------------------
      if (debug_ijt) then
        call don_cm_output_member_tends_k  &
             (nlev_hires, ncc_kou, diag_unit, Param, tcc, dpf,  &
              efchr, emfhr, cld_press, ermesg, error) 
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (error /= 0 ) return
      endif

!------------------------------------------------------------------


end subroutine don_cm_compute_vert_fluxes_k



!#####################################################################

subroutine don_cm_output_member_tends_k    &
         (nlev_hires, ncc_kou, diag_unit, Param, tcc, dpf,   &
          efchr, emfhr, cld_press, ermesg, error) 

!----------------------------------------------------------------------
!    subroutine don_cm_output_member_tends_k prints out
!    the temperature (ctfhr) and moisture (cmfhr) tendency profiles 
!    produced by the current ensemble member.
!----------------------------------------------------------------------

use donner_types_mod, only : donner_param_type

implicit none 

!----------------------------------------------------------------------
integer,                     intent(in)  :: nlev_hires, ncc_kou,  &
                                            diag_unit
type(donner_param_type),     intent(in)  :: Param
real, dimension(nlev_hires), intent(in)  :: tcc, dpf, efchr, emfhr, &
                                            cld_press
character(len=*),            intent(out) :: ermesg
integer,                     intent(out) :: error

!---------------------------------------------------------------------
!   intent(in) variables:
!
!     ncc_kou        vertical index on cloud model grid of the level
!                    above cloud top (maximum number of vertical
!                    levels affected by current cloud)
!     diag_unit      unit number for column diagnostics output, if 
!                    diagnostics are requested for the current column
!     tcc            temperature field at cloud model levels
!                    [ degrees K ]
!     dpf            condensation rate profile 
!                    [ kg(h2o) / ( kg(air) sec) ]
!     efchr          profile of entropy tendency due to flux 
!                    convergence [ deg K / sec ]
!     emfhr          profile of moisture tendency due to flux 
!                    convergence [ kg(h2o) / (kg(air) sec) ]
!     cld_press      pressure at cloud model levels [ Pa ]
!     
!----------------------------------------------------------------------

!---------------------------------------------------------------------
!   local variables:

      real, dimension(nlev_hires)     :: cmfhr, ctfhr
      real                            :: convrat
      integer                         :: k

!---------------------------------------------------------------------
!   local variables:
!
!      cmfhr     moisture tendency due to flux convergence and 
!                condensation [ kg(h2o) / (kg(air) sec) ]
!      ctfhr     entropy tendency due to flux convergence and 
!                condensation [ deg K / sec ]
!      convrat   latent heat factor used with condensation term in 
!                entropy equation [ deg K ]
!      k         do-loop index
!
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!---------------------------------------------------------------------
!    define the tendency terms at the levels within the cloud.
!---------------------------------------------------------------------
      do k=1,ncc_kou

!---------------------------------------------------------------------
!    define the appropriate latent heat to be used at the current 
!    cloud level. use the latent heat of vaporization for temperatures
!    at or above the freezing point (Param%tfre) and the latent heat of 
!    sublimation at temperatures lower than the freezing point. 
!---------------------------------------------------------------------
        if (tcc(k) >= Param%tfre) then
          convrat = Param%hlv/Param%cp_air   
        else
          convrat = Param%hls/Param%cp_air   
        endif

!---------------------------------------------------------------------
!    define the entropy tendency as the sum of the vertical flux conver-
!    gence of entropy (efchr) and the latent heat release.
!---------------------------------------------------------------------
        ctfhr(k) = -dpf(k)*convrat + efchr(k)

!---------------------------------------------------------------------
!     define the moisture tendency as the sum of vertical flux conver-
!     gence and total condensation. 
!---------------------------------------------------------------------
        cmfhr(k) = dpf(k) + emfhr(k)
      end do

!---------------------------------------------------------------------
!    set the values at layers above cloud top to 0.0.
!---------------------------------------------------------------------
      cmfhr(ncc_kou+1:nlev_hires) = 0.
      ctfhr(ncc_kou+1:nlev_hires) = 0.

!---------------------------------------------------------------------
!  ctfhr : cloud ensemble entropy tendency due to flux convergence and
!          condensation       
!          [ deg K / sec ]
!---------------------------------------------------------------------
      do k=1,nlev_hires       
        if (ctfhr(k) /= 0.0) then
          write (diag_unit, '(a, i4, 2e20.12)')  &
                    'in mulsub: k, P & cond/efc              =',    &
                     k, cld_press(k),ctfhr(k)
        endif
      end do

!---------------------------------------------------------------------
!  cmfhr : cloud ensemble moisture tendency due to flux convergence and
!          condensation       
!          [ kg(h2o) / (kg(air) sec) ]
!---------------------------------------------------------------------
      do k=1,nlev_hires       
        if (cmfhr(k) /= 0.0) then
          write (diag_unit, '(a, i4, 2e20.12)')  &
                   'in mulsub: k, P & cond/mfc              =',    &
                     k, cld_press(k),cmfhr(k)
        endif
      end do

!---------------------------------------------------------------------


end subroutine don_cm_output_member_tends_k


!#####################################################################


subroutine don_cm_process_condensate_k     &
       (nlev_lsm, nlev_hires, ntr, cldtop_indx, diag_unit, debug_ijt, &
        Param, acpre, accond, pb, pt_kou, pf, pftr, tcc, rcl, &
        cld_press, phalf_c, conint, dint, pmel, pmelt_lsm, precip, cu,&
        cell_precip, sumlhr, summel, dpf, dpftr, dfr, cell_melt, ermesg, error)

!----------------------------------------------------------------------
!
!----------------------------------------------------------------------

use donner_types_mod, only : donner_param_type

implicit none

!----------------------------------------------------------------------
integer,                      intent(in)    :: nlev_lsm, nlev_hires, &
                                               ntr, &
                                               cldtop_indx, diag_unit
logical,                      intent(in)    :: debug_ijt  
type(donner_param_type),      intent(in)    :: Param
real,                         intent(in)    :: acpre, accond, pb, pt_kou
real, dimension(nlev_hires),  intent(in)    :: pf, tcc, rcl, cld_press
real, dimension(nlev_hires,ntr), intent(in) :: pftr
real, dimension(nlev_lsm+1),  intent(in)    :: phalf_c
real,                         intent(inout) :: conint, dint,         &
                                               pmel, precip, cu, &
                                               cell_precip, sumlhr, &
                                               summel
real,                          intent(in) :: pmelt_lsm
real, dimension(nlev_hires),  intent(inout) :: dfr
real, dimension(nlev_hires),  intent(out)   :: dpf
real, dimension(nlev_hires,ntr),  intent(out)   :: dpftr
real, dimension(nlev_lsm),    intent(out)   :: cell_melt
character(len=*),             intent(out)   :: ermesg
integer,                      intent(out)   :: error
              

      real    :: dmela
      integer  :: k, kc, n
 
!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

      dpf(1:nlev_hires) = 0.
      dpftr(:,:) = 0.
      precip = 0.
      conint = 0.


      if (cldtop_indx /= 0) then
!--------------------------------------------------------------------
!    loop over the cloud levels to adjust various outputs as needed.
!--------------------------------------------------------------------
        do k=1,cldtop_indx

!--------------------------------------------------------------------
!    convert the layer-mean values of cloud-area-weighted condensation 
!    (pf) and wet deposition (pftr)
!    to values at cloud model interfaces (dpf and dpftr).
!--------------------------------------------------------------------
          if (k == 1) then
            dpf(1) =     pf(1)
            do n=1,ntr
              dpftr(1,n) =     pftr(1,n)
            end do
          else 
            dpf(k) = 0.5*(pf(k) + pf(k-1))
            do n=1,ntr
              dpftr(k,n) = 0.5*(pftr(k,n) + pftr(k-1,n))
            end do
          endif 
          if (k == cldtop_indx) then
            dpf(cldtop_indx+1) = 0.5*pf(cldtop_indx)
            do n=1,ntr
              dpftr(cldtop_indx+1,n) = 0.5*pftr(cldtop_indx,n)
            end do
          endif
        end do
      endif ! (cldtop_indx /= 0)

!--------------------------------------------------------------------
!    if in a diagnostics column, output the freezing (dfr) and conden-
!    sation (dpf) rates and cloud radius (rcl) at each cloud model
!    level, before normalization by the cloud base area.
!----------------------------------------------------------------------
      do k=1,cldtop_indx + 2
        if (debug_ijt) then
          write (diag_unit, '(a, i4, 3e20.12)')  &
                'in mulsub: k,dfr,dpr,rcl= ', k, dfr(k) ,dpf(k), rcl(k)
        endif
        if (debug_ijt) then
          if ( k > 1 ) then 
            write (diag_unit, '(a, i4, 3e20.12)')  &
                  'in mulsub: k,dpf,pf,pfavg= ', k, dpf(k) ,pf(k),  &
! note pf(k-1) when k=1 is garbage !
                                    0.5*(pf(k) + pf(k-1))
          else
            write (diag_unit, '(a, i4, 3e20.12)')  &
                  'in mulsub: k,dpf,pf,pfavg= ', k, dpf(k) ,pf(k),  &
! note pf(k-1) when k=1 is garbage !
                                    pf(k)
          endif
        endif
 
!--------------------------------------------------------------------
!    normalize the freezing rate and condensation rate profiles by the
!    cloud base area of ensemble member #1.
!--------------------------------------------------------------------
        dfr(k) = dfr(k)/(Param%cloud_base_radius**2)
        dpf(k) = dpf(k)/(Param%cloud_base_radius**2)
        do n=1,ntr
          dpftr(k,n) = dpftr(k,n)/(Param%cloud_base_radius**2)
        end do
 
!----------------------------------------------------------------------
!    if in a diagnostics column, output the freezing (dfr) and conden-
!    sation (dpf) rates and cloud radius (rcl) at each cloud model
!    level, after normalization by the cloud base area.
!----------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, i4, 3e20.12)')  &
               'in mulsub: k,dfr,dpr,rcl= ', k, dfr(k), dpf(k), rcl(k)
          do n=1,ntr
            write (diag_unit, '(a, i4, e20.12)')  &
                'in mulsub: k,dpftr= ', k, dpftr(k,n)
          end do
        endif
      end do

      if (cldtop_indx /= 0) then
        do k=1,cldtop_indx
          conint = conint + pf(k)*Param%dp_of_cloud_model/Param%grav  
        if (debug_ijt) then
          write (diag_unit, '(a, i4, 2e20.12)')  &
                          'in mulsub: kc, PF,coninT= ',k, pf(k)*  &
                                Param%dp_of_cloud_model/Param%grav, &
                                                            conint  
        endif
        end do

!--------------------------------------------------------------------
!    if in diagnostics column, output the ensemble member number (kou), 
!    moisture convergence integral (sub1sum),  and condensation 
!    integral (conint).
!--------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, e20.12, a)')  &
                   'in cloudm: CONINT= ',CONINT,' KG/(M**2)/SEC'
        endif

!--------------------------------------------------------------------
!    normalize the column integral of condensate by the cloud base area
!    of ensemble member #1. 
!--------------------------------------------------------------------
        conint = conint/(Param%cloud_base_radius**2)

!--------------------------------------------------------------------
!    define the precipitation generated by this ensemble member.
!--------------------------------------------------------------------
        precip = conint*acpre/accond
      endif ! (cldtop_indx /= 0)

!--------------------------------------------------------------------
!    define the condensation (cu) and precipitation rates (cell_precip) 
!    in units of mm(h2o) per day. add this ensemble member's contribution
!    to the total precipitation (ensmbl_precip) and condensation
!    (ensmbl_cond), normalized by the ensemble member's cloud base area.
!--------------------------------------------------------------------
      cu  = conint*Param%seconds_per_day
      cell_precip  = precip*Param%seconds_per_day
 
!----------------------------------------------------------------------
!    if in a diagnostics column, output the condensation (conint),
!    precipitataion (precip) and freezing (dint) rates for this ensemble
!    member, in both kg per m**2 per second and in mm per day.
!----------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 3e20.12, a)')  &
                   'in mulsub: CONPRE, CPRE,DINT= ', conint, precip, &
                      dint, ' KG/(M**2)/SEC'
        write (diag_unit, '(a, e20.12, a)')   &
                    'in mulsub: CONDENSATION PRE= ', cu, ' MM/DAY'
        write (diag_unit, '(a, e20.12, a)')  &
                 'in mulsub: CLOUD MODEL PRE= ', cell_precip, ' MM/DAY'
      endif

!---------------------------------------------------------------------
!    compute the contribution to the column integrals of total conden-
!    sation (sumlhr) and frozen condensate from the lowest cloud layer.
!    units are kg(h2o) per square meter per second.
!---------------------------------------------------------------------
      sumlhr = (dpf(1)*((Param%dp_of_cloud_model/2.)/Param%grav))
      if (tcc(1) <= Param%tfre) then
        summel = (Param%dp_of_cloud_model/2.)*dpf(1)/Param%grav  
      else
        summel = 0.
      endif

!---------------------------------------------------------------------
!    if in a diagnostics window, output the condensation rate and
!    accumulated condensate integral at level 1.
!---------------------------------------------------------------------
      if (debug_ijt) then
        kc = 1
        write (diag_unit, '(a, i4, 2e20.12)')  &
                               'in mulsub: kc,DPF, SUMLHR= ',kc, &
               dpf(1)*((Param%dp_of_cloud_model/2.)/Param%grav),   &
                                                sumlhr
      endif

!--------------------------------------------------------------------
!    add the contributions from the remaining cloud layers to the 
!    integrals.
!--------------------------------------------------------------------
      do kc = 2,cldtop_indx+1

!---------------------------------------------------------------------
!    add the density weighted increment of column condensate to the 
!    array accumulating it (sumlhr). if in diagnostic column, write
!    out the condensate at the level and the current vertical sum of 
!    the column condensate.
!---------------------------------------------------------------------
        sumlhr = sumlhr + (dpf(kc)*(Param%dp_of_cloud_model/Param%grav))
        if (debug_ijt) then
          write (diag_unit, '(a, i4, 2e20.12)')  &
                         'in mulsub: kc,DPF,SUMLHR= ',kc,   &
               dpf(kc)*(Param%dp_of_cloud_model/Param%grav),   &
                                                               sumlhr
        endif

!---------------------------------------------------------------------
!    if the temperature is at or below the freezing level, add the 
!    current-level's density-weighted condensation to the array accum-
!    ulating the total frozen condensate in the column (summel).
!---------------------------------------------------------------------
        if (tcc(kc) <= Param%tfre) then
          summel = summel +  Param%dp_of_cloud_model*dpf(kc) /Param%grav 
        endif
      end do

!---------------------------------------------------------------------
!    if in diagnostics column, output the column integral of frozen 
!    condensate (summel). determine the amount of this condensate which
!    precipitates out by multiplying by the factor (cell_precip/cu).  add
!    this amount to the amount frozen in the updraft (dint) to obtain the
!    total amount which must be melted in the column (dints).
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12)')  'summel= ',summel
      endif

!--------------------------------------------------------------------
!    define the cloud model pressure at the level where melting occurs 
!    (pmel). 
!---------------------------------------------------------------------
      pmel = pb             
      do kc=1,cldtop_indx
        if ((tcc(kc) >= Param%kelvin) .and.   &
            (tcc(kc+1) <= Param%kelvin))  pmel = cld_press(kc)
      end do
      if (tcc(cldtop_indx+1) == Param%kelvin) pmel = pt_kou 
          if (debug_ijt) then
            write (diag_unit, '(a, 3e20.12 )')  &
             'pmelt_lsm, pmel from hi-res, pb   = ', pmelt_lsm, pmel, pb
          endif

!---------------------------------------------------------------------
!    if there has been no freezing in the column, then there will be
!    no melting.
!---------------------------------------------------------------------
      dmela = 0.
      if (dint == 0.) then
        cell_melt(1:nlev_lsm) = 0.
      else

!--------------------------------------------------------------------
!    if the melting level is above cloud base, partition the integrated
!    ice melt from the cloud model within the appropriate large-scale 
!    model layers.
!--------------------------------------------------------------------
        if (pb > pmelt_lsm) then

!--------------------------------------------------------------------
!    define the rate of ice melt (dmela) in units of g(h2o) per 
!    kg(air) per day. 
!--------------------------------------------------------------------
          if (cu /= 0.0) then
! melt ice precip plus frozen liq precip
            dmela = -((summel*cell_precip/cu + dint*cell_precip/cu)*  &
                                  Param%grav/(pmelt_lsm - pb))*8.64E07
         endif

!--------------------------------------------------------------------
!    call map_hi_res_intgl_to_lo_res_col to distribute this integrated 
!    column melting over the appropriate pressure interval on the 
!    large-scale model grid. if in a diagnostic column, output the 
!    melting rates as mapped to the large-scale grid.
!--------------------------------------------------------------------
          if (debug_ijt) then
            write (diag_unit, '(a, e20.12, 2f19.10)')  &
              'in cm_intgl_to_gcm_col: dmela,pb,pmelt_lsm= ',  &
                                                 dmela, pb, pmelt_lsm
          endif
          call don_u_map_hires_i_to_lores_c_k  &
            (nlev_lsm, dmela, pb, pmelt_lsm, phalf_c, cell_melt, ermesg, error)
 
!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
          if (error /= 0 ) return

!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
          if (debug_ijt) then
            do k=1,nlev_lsm              
              if (cell_melt(k) /= 0.0) then
                write (diag_unit, '(a, i4, e20.12)') &
                 'in cm_intgl_to_gcm_col: k,cell_melt= ',k,cell_melt(k)
              endif
            end do
          endif

!---------------------------------------------------------------------
!    if the melting level is at or below cloud base, then there is no
!    in-cloud melting of frozen condensate.
!---------------------------------------------------------------------
        else
          cell_melt(1:nlev_lsm) = 0.
        endif  ! ( pb > pmel)
      endif  !  (dint == 0.0)


!---------------------------------------------------------------------
!    change the sign of the integral and convert to units of kg(h2o)
!    per square meter per day (or mm per day). if in diagnostics col-
!    umn, output the column integrated melting (summel), precipitation 
!    (cell_precip) and condensation rates (cu) in units of mm per day.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12,a)')   &
             'in mulsub: rc,cu= ',    &
        cell_precip,cu, &
                                            'mm/day'
        if ( cu .ne. 0.0 ) then
          write (diag_unit, '(a, e20.12,a)')   &
             'in mulsub: summel = ',    &
        -(summel*cell_precip/cu)*Param%seconds_per_day, 'mm/day'
        else
          write (diag_unit, '(a)')   &
             'in mulsub: summel = undefined due to cu = 0'
        endif

      endif

!---------------------------------------------------------------------


end subroutine don_cm_process_condensate_k



!#####################################################################


subroutine don_cm_simult_k   &
         (diag_unit, debug_ijt, lfc_not_reached, Param, pcsave, rmu,  &
          cloud_temp, cloud_radius, w_vel, cloud_p, liq_wat, env_temp, &
          env_mixing_ratio, dtdp,  drdp, dwdp, ermesg, error)

!--------------------------------------------------------------------
!    subroutine simult returns the vertical derivatives of temperature
!    (dtdp), cloud radius (drdp) and vertical velocity (dwdp) in the
!    cloud.
!    Reference: Donner, JAS, 1986, v43, pp.2277-2288.
!    See LJD "Cloud Model 89" notes on Generalized mu (10/1/89)
!    and dwdz (10/2/89). The value of epsilon is taken as 1.
!    Version where cloud properties independent of cloud area.
!--------------------------------------------------------------------

use donner_types_mod, only : donner_param_type
use sat_vapor_pres_k_mod, only: compute_qs_k

implicit none

!--------------------------------------------------------------------
integer,                 intent(in)    ::  diag_unit
logical,                 intent(in)    ::  debug_ijt, lfc_not_reached
type(donner_param_type), intent(in)    ::  Param
real,                    intent(in)    ::  pcsave, rmu, cloud_temp,  &
                                           cloud_radius, w_vel,   &
                                           cloud_p, liq_wat, env_temp, &
                                           env_mixing_ratio
real,                    intent(out)   ::  dtdp, drdp, dwdp
character(len=*),        intent(out)   ::  ermesg
integer,                 intent(out)   ::  error

!---------------------------------------------------------------------
!   intent(in) variables:
!
!        cloud_temp  cloud temperature [ deg K ]
!        cloud_radius 
!                    cloud radius  [ m ]
!        w_vel       cloud vertical velocity [ m / s ] 
!        cloud_p     pressure at current level [ Pa ]
!        env_temp    environmental temperature [ deg K ]
!        env_mixing_ratio
!                    environmental mixing ratio [ kg / kg ]
!        liq_wat     liquid water [ kg / kg ]
!        pcsave      pressure at which cloud ensemble member 1 becomes
!                    buoyant [ Pa ]
!        rmu         entrainment coefficient [ m**(-1) ]
!        lfc_not_reached   if true, have not yet reached buoyant part of 
!                    cloud (vertical velocity has not yet begun to
!                    increase)
!        debug_ijt   is this a diagnostics column ?
!        diag_unit   output unit number for this diagnostics column
!
!   intent(out) variables:
!
!        dtdp        temperature derivative [ deg K / Pa ]
!        drdp        cloud radius derivative [ m / Pa ]
!        dwdp        vertical velocity derivative [ m / (sec Pa) ]
!
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!   local variables:

      real  ::  areap, c2, c3, c4, dadp, dw2dz, dwdz, dzdp, dz, es,  &
                entrain_co, htv, htve, lat, rhoaw, rhoawp, &
                rstar, sphum, tcae, teae, test, west, w2test, wtest
      integer :: nbad

!---------------------------------------------------------------------
!   local variables:
!
!         areap             updated cloud area [ m**2 ]
!         c2                term in dt/dz eqn (eqn 5, Donner, JAS, 1993)
!         c3                term in dt/dz eqn (eqn 5, Donner, JAS, 1993)
!         c4                term in dt/dz eqn (eqn 5, Donner, JAS, 1993)
!         dadp              derivative of cloud area with respect to
!                           pressure [ m**2 / Pa ]
!         dw2dz             vertical derivative of vertical velocity
!                           squared with respect to height 
!                           [ m / sec**2 ]
!         dwdz              vertical derivative of vertical velocity
!                           with respect to height [ sec**(-1) ]
!         dzdp              dz/dp when virtual mass coefficient is
!                           used [ m / Pa ]
!         dz                delta z corresponding to dp_of_cloud_model
!                           [ m ]
!         es                saturation vapor pressure [ Pa ]
!         entrain_co        entrainment coefficient [ m**(-1) ]
!         htv               virtual temperature of cloud       [ deg K ]
!         htve              virtual temperature of environment [ deg K ]
!         lat               latent heat relevant at input temperature
!                           cloud_temp [ J/kg ]
!         rhoaw             rho * cloudarea * vertical velocity for
!                           input values [ kg / sec ]
!         rhoawp            rho * cloudarea * vertical velocity for
!                           updated values [ kg / sec ]
!         rstar             gas constant for moist air [ J / (kg degK) ]
!         sphum             saturation specific humidity [ kg / kg ]   
!         tcae              equivalent potential temperature in cloud
!                           [ deg K ]
!         teae              equivalent potential temperature in env-
!                           ironment [ deg K ]
!         test              updated temperature [ deg K ]
!         west              updated vetical velocity [ m / sec ]
!         w2test            updated vertical velocity squared (without
!                           entrainment term) [ m**2 / sec**2 ]
!         wtest             updated vertical velocity  (without
!                           entrainment term) [ m / sec ]
!
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!---------------------------------------------------------------------
!   define appropriate latent heat for the given cloud temperature.
!---------------------------------------------------------------------
      if (cloud_temp < Param%tfre) then
        lat = Param%hls    
      else
        lat = Param%hlv   
      endif
 
!---------------------------------------------------------------------
!    define the specific humidity at the cloud temperature (sphum) and
!    the gas constant for moist air (rstar).
!---------------------------------------------------------------------
      call compute_qs_k (cloud_temp, cloud_p, Param%D622,  &
                           Param%D608, sphum, nbad, esat = es)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (nbad /= 0) then
        ermesg = 'subroutine don_cm_simult_k: '// &
                 'temperatures out of range of esat table'
        error = 1
        return
      endif

      rstar = Param%rdgas*(1. + Param%d608*sphum)

!---------------------------------------------------------------------
!    define the in-cloud (htv) and environmental (htve) virtual 
!    temperature.                                                       
!---------------------------------------------------------------------
      htv  = cloud_temp*(1. + Param%d608*sphum)
      htve = env_temp*(1. + Param%d608*(env_mixing_ratio/   &
                                                (1.0+env_mixing_ratio)))

!---------------------------------------------------------------------
!    if in diagnostics window, output the in-cloud and environmmental
!    virtual temperatures.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)')  &
                      'in simult: htve,htv= ',htve, htv
      endif

!---------------------------------------------------------------------
!    define one of the terms present in the dt/dz eqn (Eqn 5) in 
!    Donner, JAS, 1993, v50, pp.890-906.
!---------------------------------------------------------------------
      c2 = (htv + (lat*sphum/rstar))*Param%grav/(Param%cp_air*cloud_temp)

!---------------------------------------------------------------------
!    if in diagnostics window, output the c2 term in the dt/dz equation
!    and the environmmental temperature, moisture and pressure at this
!    level. 
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12)')  'in simult: c2= ',c2
        write (diag_unit, '(a, f20.14, f19.10, e20.12)')  &
              'in simult: te,p,qe= ',env_temp, cloud_p ,env_mixing_ratio
      endif

!--------------------------------------------------------------------
!    call tae to calculate environmental adiabatic equivalent 
!    temperature (teae).
!--------------------------------------------------------------------
      call don_cm_tae_k    &
           (Param, env_temp, cloud_p, env_mixing_ratio, lat, teae,  &
            ermesg, error)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (error /= 0 ) return

!---------------------------------------------------------------------
!    if in diagnostics window, output the environmental adiabatic 
!    equivalent temperature.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, f20.14)') 'in simult: teae= ',teae
      endif

!--------------------------------------------------------------------
!    calculate in-cloud adiabatic equivalent temperature (tcae).
!--------------------------------------------------------------------
      tcae = cloud_temp*exp(lat*sphum/(Param%cp_air*cloud_temp))

!----------------------------------------------------------------------
!    define dz/dp in the presence of a virtual mass coefficient (used
!    to roughly account for non-hydrostatic effects). define the cor-
!    resonding dz for the cloud model pressure increment dp.
!    Reference: Donner, JAS, 1986, v43, pp.2277-2288.
!----------------------------------------------------------------------
      dzdp = -Param%rdgas*(1. + Param%virt_mass_co)*htv*htve/    &
                                (Param%grav*cloud_p*  &
                                      (Param%virt_mass_co*htve + htv))
      dz   = Param%dp_of_cloud_model*dzdp

!---------------------------------------------------------------------
!    define the remaining terms (c3, c4) in the dt/dz equation (eqn (5)
!    in the reference cited below).
!    note : the entrainment coefficient (dln(rho*a*w)/dz) is given by 
!    (exp(rmu*dz) - 1.)/(exp(rmu*dz)*dz), with assumption that 
!    rhoaw(z) = rhoaw(0)*exp(rmu*dz)
!    Reference: Donner, JAS, 1993, v50, pp.890-906.
!---------------------------------------------------------------------
      entrain_co = (exp(rmu*dz) - 1.0)/(exp(rmu*dz)*dz)
      c3 = (tcae - teae)/exp(lat*sphum/(Param%cp_air*cloud_temp))*  &
                                                               entrain_co
      c4 = 1. + (Param%d622*lat*(Param%hlv*es/     &
                   (Param%rvgas*(cloud_temp**2)))/(Param%cp_air*cloud_p))

!---------------------------------------------------------------------
!    if in diagnostics window, output the c3 and c4 terms in the dt/dz 
!    equation.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)')  'in simult: c3,c4= ',c3,c4
      endif

!---------------------------------------------------------------------
!    define the dt/dp tendency by combining the component terms 
!    in the dt/dz equation and multiplying by dz/dp. define an 
!    updated estimated value of in-cloud temperature (test).
!---------------------------------------------------------------------
      dtdp = -((c2 + c3)/c4)*dzdp
      test = cloud_temp + dtdp*Param%dp_of_cloud_model

!--------------------------------------------------------------------
!    define the buoyancy and drag terms in the d(w**2)/dz equation. 
!    Reference: Donner, JAS, 1993, v50, pp.890-906.
!--------------------------------------------------------------------
      dw2dz = 2.0*((Param%grav*(htv - htve)/      &
                                  (htve*(1. + Param%virt_mass_co))) -  &
                                                      Param%grav*liq_wat)

!-------------------------------------------------------------------
!    produce an updated vertical velocity squared (wtest) by adding this
!    tendency term to the input value. be sure w**2 is positive; then 
!    define the vertical velocity.
!-------------------------------------------------------------------
      w2test = (w_vel**2) + dw2dz*dz
      if (w2test < 0.) then
        wtest = 0.
      else
        wtest = sqrt(w2test)
      endif

!---------------------------------------------------------------------
!    if in diagnostics window, output the updated vertical velocity.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, l4, e20.12)')  &
                     'in simult: testlc,wtest= ',lfc_not_reached, wtest
      endif

!----------------------------------------------------------------------
!    define the dw/dp derivative from this updated value, ignoring the 
!    entrainment term.
!----------------------------------------------------------------------
      dwdp = (wtest - w_vel)/Param%dp_of_cloud_model

!----------------------------------------------------------------------
!    determine if the necessary conditions for entrainment are present.
!    allow entrainment to change vertical velocity if parcel has 
!       (1) previously achieved initial acceleration 
!           (.not. lfc_not_reached) 
!           and it has ascended above the level where the most entrain-
!           ing parcel (ensemble member #1) initially accelerates or,
!       (2) the parcel is currently accelerating. 
!    the first condition ensures that a shallow, lower-entrainment cloud
!    will not develop below the pressure where the most entraining
!    parcel initially accelerates.
!----------------------------------------------------------------------
      if ( ((.not. lfc_not_reached) .and. (cloud_p <= pcsave)) .or.   &
           (dwdp <= 0.) ) then
        dwdz = -w_vel*entrain_co
        dwdp = (dwdz*dzdp) + dwdp
      endif

!----------------------------------------------------------------------
!    if the parcel has not reached the level of free convection, do
!    not allow its vertical velocity to decrease. BL turbulence or 
!    other sub-grid mechanisms are assumed present to sustain its
!    upward motion.
!----------------------------------------------------------------------
      if ((lfc_not_reached) .and. (dwdp > 0.) )  then
        dwdp = 0.
      endif

!----------------------------------------------------------------------
!    if the parcel is above its level of free convection  but below
!    the pressure level at which the ensemble member #1 initially 
!    accelerates, do not allow its vertical velocity to decrease. this
!    ensures that less entraining clouds will not develop between the
!    ground and the pressure at which the most entraining parcel 
!    initally accelerates. BL turbulence or some other sub-grid 
!    mechanism is assumed to sustain its upward motion.
!----------------------------------------------------------------------
      if ( (.not. lfc_not_reached) .and. (cloud_p > pcsave) .and.    &
           (dwdp > 0.) )  then
        dwdp = 0.
      endif

!---------------------------------------------------------------------
!    if in diagnostics window, output the updated vertical velocity.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a,l4,  2e20.12)')  &
          'in simult: testlc,dwdp,test= ',lfc_not_reached,dwdp,test
      endif

!---------------------------------------------------------------------
!    produce an updated vertical velocity (eqn (6) of reference).
!---------------------------------------------------------------------
      west = w_vel + dwdp*Param%dp_of_cloud_model

!--------------------------------------------------------------------
!    calculate the terms in the cloud area tendency equation.
!    Reference: Donner, JAS, 1993, v50, pp.890-906.
!--------------------------------------------------------------------
      if (west < Param%wdet) then

!--------------------------------------------------------------------
!    if the updated vertical velocity is weaker than the detrainment
!    velocity, set the cloud radius tendency to 0.0,maintaining the
!    input value.
!--------------------------------------------------------------------
        drdp = 0.
      else

!--------------------------------------------------------------------
!    define the value of rho*cloudarea*vertical velocity (rhoaw) for 
!    the input values of these variables. define the same quantity at 
!    a level dz higher (rhoawp), assuming rhoaw(z) = 
!    rhoaw(0)*exp(rmu*dz).
!--------------------------------------------------------------------
        rhoaw = cloud_p*(cloud_radius**2)*w_vel/(Param%rdgas*htv)
        rhoawp = rhoaw*exp((dzdp*rmu)*Param%dp_of_cloud_model)

!---------------------------------------------------------------------
!    if in diagnostics window, output rho*cloudarea*w obtained from
!    input values and at a height dz higher.
!---------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, f19.10, 2e20.12)')  &
                       'in simult: p,fm,fmp= ', cloud_p, rhoaw, rhoawp
        endif

!---------------------------------------------------------------------
!    define the virtual temperature (htv) corresponding to the updated 
!    temperature test.
!----------------------------------------------------------------------
        call compute_qs_k (test, cloud_p, Param%D622,  &
                           Param%D608, sphum, nbad)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (nbad /= 0) then
          ermesg = 'subroutine don_cm_simult_k: '// &
                   'temperatures out of range of esat table'
          error = 1
          return
        endif
        htv = test*(1. + Param%d608*sphum)

!----------------------------------------------------------------------
!    define the updated cloud area (areap). define the area tendency 
!    (dadp) and then the cloud radius tendency (drdp).
!----------------------------------------------------------------------
        areap = rhoawp*Param%rdgas*htv/  &
                               ((cloud_p + Param%dp_of_cloud_model)*west)
        dadp = (areap - (cloud_radius**2))/Param%dp_of_cloud_model
        drdp = dadp/(2.*cloud_radius)
      endif

!--------------------------------------------------------------------



end subroutine don_cm_simult_k




!####################################################################

subroutine don_cm_tae_k    &
         (Param, init_temp, init_pr, parcel_mixing_ratio, latent_heat,  &
          equivalent_temp, ermesg, error)

!--------------------------------------------------------------------
!    subroutine tae determines the saturation temperature of an init-
!    ially non-saturated parcel ( defined by init_temp, 
!    parcel_mixing_ratio, init_pr) by incrementally moving the parcel 
!    upward dry adiabatically until it reaches the temperature at which 
!    saturation would occur for the initial pressure init_pr.  using this
!    saturation temperature (te) and the parcel moisture 
!    parcel_mixing_ratio and pressure init_pr, the equivalent temperature
!    is calculated (the temperature attained if all vapor were condensed 
!    out by moist adiabatic ascent, followed by dry adiabatic descent to
!    the initial pressure level init_pr.
!--------------------------------------------------------------------

use donner_types_mod, only : donner_param_type
use sat_vapor_pres_k_mod, only: compute_mrs_k

implicit none

!--------------------------------------------------------------------
type(donner_param_type), intent(in)    :: Param
real,                    intent(in)    :: init_temp, init_pr,   &
                                          parcel_mixing_ratio,  &
                                          latent_heat
real,                    intent(out)   :: equivalent_temp
character(len=*),        intent(out)   :: ermesg
integer,                 intent(out)   :: error

!--------------------------------------------------------------------
!   intent(in) variables:
!
!        init_temp              temperature   [ deg K ]
!        init_pr                pressure      [ Pa ]
!        parcel_mixing_ratio    mixing ratio [kg(H2O)/kg (dry air) ]
!        latent_heat            applicable latent heat constant  [ J/kg ]
!
!   intent(out) variables:
!
!        equivalent_temp        adiabatic equivalent temperature 
!                               [ deg K ]
!
!---------------------------------------------------------------------

!--------------------------------------------------------------------
!   local variables:

      real      ::   pr, te, mre
      integer   ::   nbad

!--------------------------------------------------------------------
!   local variables:
!
!        pr   pressure at current level  [ Pa ]
!        te   temperature of parcel move dry adiabatically from
!             pressure p to pressure pr
!        es   stauration vapor pressure at temperature te
!        mre  saturation mixing ratio at temperature te and pressure p
!        k    do-loop index
!
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutinE.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!---------------------------------------------------------------------
!    initialize the output variable equivalent_temp and the pressure at 
!    current level (pr).
!---------------------------------------------------------------------
      equivalent_temp = init_temp
      pr = init_pr

!---------------------------------------------------------------------
!    vertical loop. if the pressure  pr is less than pstop (absolute 
!    lowest pressure for cloud) exit the loop.
!---------------------------------------------------------------------
      do while (pr >= Param%pstop) 

!--------------------------------------------------------------------
!    define the temperature (te) at pressure pr assuming dry adiabatic
!    ascent from initial pressure p. determine the saturation vapor
!    pressure (es) for this temperature. define the saturation mixing
!    ratio for a parcel of this temperature at the initial pressure
!    level (mre).
!--------------------------------------------------------------------
        te = init_temp*((pr/init_pr)**Param%kappa)
        call compute_mrs_k (te, init_pr,                               &
                            Param%d622 , Param%d608 , mre, nbad, &
                            mr = parcel_mixing_ratio)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
        if (nbad /= 0) then
          ermesg = 'subroutine don_cm_tae_k: '// &
                   'temperatures out of range of esat table'
          error = 1
          return
        endif

!--------------------------------------------------------------------
!    determine if saturation at the initial pressure level would occur
!    for this temperature.
!--------------------------------------------------------------------
        if (parcel_mixing_ratio >= mre) then

!--------------------------------------------------------------------
!    if saturation would occur for this temperature (te) at the initial
!    pressure init_pr (i.e., parcel_mixing_ratio >= mre), then use te to
!    calculate the equivalent temperature. exit the loop. otherwise,
!    increment the current pressure and continue within the loop.
!--------------------------------------------------------------------
          equivalent_temp = init_temp*exp(latent_heat*  &
                                   parcel_mixing_ratio/(te*Param%cp_air))
          exit
        endif
        pr = pr + Param%dp_of_cloud_model
      end do

!---------------------------------------------------------------------



end subroutine don_cm_tae_k

!####################################################################

subroutine don_cm_micro_k   &
         (diag_unit, debug_ijt, Param, tc1, tc2, p1, p2, te1, te2,  &
          qe1, qe2, w1, w2, rr, rmu, qrw, qcw, qlw, dcw1, dqrw3, ermesg, error)

!----------------------------------------------------------------------
!    subroutine micro calculates microphysical tendencies (kessler
!    microphysics) of the cloud parcel defined by (tc, w, qrw, qcw, 
!    qlw, rr) in an environment defined by (te, qe, rmu) during the 
!    movement of the parcel from pressure level p1 to level p2.
!--------------------------------------------------------------------

use donner_types_mod, only : donner_param_type
use sat_vapor_pres_k_mod, only: compute_mrs_k

implicit none

!--------------------------------------------------------------------
integer,                  intent(in)    :: diag_unit
logical,                  intent(in)    :: debug_ijt
type(donner_param_type),  intent(in)    :: Param
real,                     intent(in)    :: tc1, tc2, p1, p2, te1, te2, &
                                           qe1, qe2, w1, w2, rr, rmu
real,                     intent(inout) :: qrw, qcw, qlw, dcw1, dqrw3
character(len=*),         intent(out)   :: ermesg
integer,                  intent(out)   :: error

!--------------------------------------------------------------------
!  intent(in) variables:
!
!        tc1       cloud temperature at starting point [ deg K ]
!        tc2       cloud temperature at ending point [ deg K ]
!        p1        pressure at starting point [ Pa ]
!        p2        pressure at ending point [ Pa ]
!        te1       environmental temperature at starting point [ deg K ]
!        te2       environmental temperature at ending point [ deg K ]
!        qe1       environmental mixing ratio at starting point  
!                  [ kg(H2O) / kg(air) ] 
!        qe2       environmental mixing ratio at ending point  
!                  [ kg(H2O) / kg(air) ] 
!        w1        cloud vertical velocity at starting point [ m / sec ]
!        w2        cloud vertical velocity at ending point [ m / sec ]
!        rr        cloud radius [ m ]
!        rmu       entrainment coefficient [ m**(-1) ]
!        debug_ijt is this a diagnostics column ?
!        diag_unit output unit number for this diagnostics column
!
!   intent(inout) variables:
!
!        qrw    rain water              [ g(h2o) / m**3 ]
!        qcw    cloud water             [ g(h2o) / m**3 ]
!        qlw    total liquid water      [ kg(h2o) / kg(air) ]
!        dcw1   condensation increment  [ g(h2o) / m**3 ]
!        dqrw3  precipitation increment [ g(h2o) / m**3 ]
!
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!   local variables:
  
      real :: cond, d1, d2, dcw2, dqcw3, dt_micro, dz, ent, pav,  &
              rb, qcwa, qeb, qrwa, red, rho, rs1, rs2, tcb, wav
      integer :: nbad

!--------------------------------------------------------------------
!   local variables:
!
!     cond      condensation in the layer [ kg(h2o) / kg(air) ]
!     d1        dz/dp at the starting level [ m / Pa ]  
!     d2        dz/dp at the ending level [ m / Pa ]  
!     dcw2      change in cloudwater due to autoconversion 
!               [ g (h2o) / m**3 ]
!     dqcw3     change in cloudwater content due to accretion     
!               [ g (h2o) / m**3 ]
!     dt_micro  microphysical time step [ sec ]
!     dz        average delta z in the layer [ m ]
!     ent       ratio of parcel mass after entrainment to mass
!               before entrainment [ dimensionless ]
!     es1       saturation vapor pressure at level p1 [ Pa ]
!     es2       saturation vapor pressure at level p2 [ Pa ]
!     pav       average pressure in the layer [ Pa ]
!     rb        average in-cloud mixing ratio in the layer 
!               [ kg(h2o) / kg(air) ]
!     qcwa      initial value of cloudwater after microphysical term
!               updates; it is then adjusted to avoid negative values
!               [ g(h2o) / m**3 ]
!     qeb       average environmental mixing ratio in the layer 
!               [ kg(h2o) / kg(air) ]
!     qrwa      initial value of rainwater after microphysical term
!               updates; it is then adjusted to avoid negative values
!               [ g(h2o) / m**3 ]
!     red       reducing factor applied to the microphysical loss terms
!               in the cloud and rain equations to avoid creation of
!               negative cloud and rain water [ dimensionless ]
!     rho       atmospheric density [ kg / m**3 ]
!     rs1       saturation mixing ratio in cloud at level p1
!     rs2       saturation mixing ratio in cloud at level p2
!     tcb       mean in-cloud temperature in the layer [ deg K ]
!     teb       mean environmental temperature in the layer [ deg K ]
!     wav       mean vertical velocity in the layer [ m / sec ]
!
!--------------------------------------------------------------------

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!--------------------------------------------------------------------
!    define the saturation mixing ratio (rs1) at level p1.
!--------------------------------------------------------------------
      call compute_mrs_k (tc1, p1, Param%d622 , Param%d608 , rs1, nbad)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (nbad /= 0) then
        ermesg = 'subroutine don_cm_micro_k: '// &
                 'temperatures out of range of esat table'
        error = 1
        return
      endif

!--------------------------------------------------------------------
!    define the saturation mixing ratio (rs2) at level p2.
!--------------------------------------------------------------------
      call compute_mrs_k (tc2, p2, Param%d622 , Param%d608 , rs2, nbad)

!----------------------------------------------------------------------
!    determine if an error message was returned from the kernel routine.
!    if so, return to calling program where it will be processed.
!----------------------------------------------------------------------
      if (nbad /= 0) then
        ermesg = 'subroutine don_cm_micro_k: '// &
                 'temperatures out of range of esat table'
        error = 1
        return
      endif

!--------------------------------------------------------------------
!    if in diagnostic column, output the relevant atmospheric fields 
!    at the two pressure levels.
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)')   &
                         'in micro: qrw,qcw= ',qrw,qcw
        write (diag_unit, '(a, e20.12)')  'in micro: rr= ',rr
        write (diag_unit, '(a, 2e20.12, 2f19.10)') &
                         'in micro: rs1,rs2,p1,p2= ',rs1,rs2,p1,p2
        write (diag_unit, '(a, 2e20.12, 2f20.14)')  &
!                        'in micro: es1,es2,tc1,tc2= ',es1,es2,tc1,tc2
                         'in micro: tc1,tc2= ',tc1,tc2
        write (diag_unit, '(a, 2f20.14)')  &
                          'in micro: te1,te2= ',te1,te2
      endif

!--------------------------------------------------------------------
!    define the layer-mean cloud temperature (tcb) and mixing ratio 
!    (rb), environmental mixing ratio (qeb), vertical velocity (wav), 
!    pressure (pav) and density (rho). define a layer-mean dz, as the  
!    average of the values of dz/dp (when using a virtual mass coeffic-
!    ient) at the two input levels. define the ratio of post-entrainment
!    mass to the pre-entrainment mass (ent).
!--------------------------------------------------------------------
      tcb = 0.5*(tc1 + tc2)
      rb  = 0.5*(rs1 + rs2)
      qeb = 0.5*(qe1 + qe2)
      wav = 0.5*(w1  +  w2)
      pav = 0.5*(p1  +  p2)
      rho = pav/(Param%rdgas*tcb*(1. + Param%d608*(rb/(1.0+rb))))
      d1  = Param%rdgas*(1. + Param%virt_mass_co)*tc1*te1/  &
            (Param%grav*p1*(Param%virt_mass_co*te1 + tc1))
      d2  = Param%rdgas*(1. + Param%virt_mass_co)*tc2*te2/  &
            (Param%grav *p2*(Param%virt_mass_co*te2 + tc2))
      dz  = -0.5*(d1 + d2)*Param%dp_of_cloud_model
      ent = exp(rmu*dz)

!--------------------------------------------------------------------
!    if in diagnostic column, output the entrainment coefficient (rmu)
!    and the environmental mixing ratio (qeb).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)')  &
                   'in micro: qeb,rmu= ', qeb, rmu
      endif

!----------------------------------------------------------------------
!    calculate the condensation (cond). the condensation is made up of
!    that which would occur due to adiabatic expansion (rs2 - rs1)
!    less the vapor which must be evaporated to saturate the environ-
!    mental air entrained into the cloud ((ent-1.)*(rs2-qeb)).  cond is
!    further modified to allow the condensation to be added to the
!    parcel at any point during the timestep through the variable
!    tr_insert_time. if the condensate is assumed to be added to the 
!    parcel while at mass m (its initial size), tr_insert_time = 0; 
!    if assumed to be added at end of step when parcel mass is ent 
!    times larger than its initial value, tr_insert_time is 1.  force 
!    the condensation to be non-negative.
!----------------------------------------------------------------------
      cond = (rs2 - rs1) + (ent - 1.)*(rs2 - qeb)
      cond = -cond/(1. - Param%tr_insert_time +    &
                         Param%tr_insert_time*ent)
      cond  = max (cond, 0.0)

!--------------------------------------------------------------------
!    if in diagnostic column, output the condensation (cond).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, e20.12)') 'in micro: cond= ', cond
      endif

!--------------------------------------------------------------------
!    convert the condensation from units of kg(h2o)/ kg(air) to units
!    of (g(h2o) / m**3) by multiplying by rho*1000. define the timestep
!    for microphysical processes (dt_micro).
!--------------------------------------------------------------------
      dcw1 = cond*rho*1000.
      dt_micro = dz/wav

!--------------------------------------------------------------------
!    calculate the amount of cloud autoconverted to rain (dcw2). 
!--------------------------------------------------------------------
      if (qcw >= Param%autoconv_threshold)  then
        dcw2 = Param%autoconv_rate*(qcw - Param%autoconv_threshold)* &
                                                                dt_micro 
      else
        dcw2 = 0.0
      endif

!--------------------------------------------------------------------
!    calculate the cloud accretion by rainwater (dqcw3).
!--------------------------------------------------------------------
      if (qcw /= 0.0 .and. qrw /= 0.0) then
        dqcw3 = 5.26e-03*qcw*(qrw**.875)*dt_micro 
      else
        dqcw3 = 0.
      endif

!--------------------------------------------------------------------
!    calculate effect of entrainment on cloud water.
!--------------------------------------------------------------------
      qcw = qcw/ent

!---------------------------------------------------------------------
!    add the microphysical terms to the cloud water equation. if neces-
!    sary adjust the magnitudes of the loss terms so that negative
!    values of cloud water are not created. define the final value of
!    qcw.
!---------------------------------------------------------------------
      qcwa = qcw + dcw1 - dcw2 - dqcw3
      if (qcwa < 0.) then
        red = (qcw + dcw1)/(dcw2 + dqcw3)
        dcw2 = dcw2*red
        dqcw3 = dqcw3*red
      endif
      qcw = qcw + dcw1 - dcw2 - dqcw3

!--------------------------------------------------------------------
!    if in diagnostic column, output the cloud water entrainment  
!    factor (ent), the updated cloud water (qcw) and the condensation 
!    term (dcw1).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 3e20.12)')  &
                             'in micro: ent,qcw,dcw1= ',ent,qcw,dcw1
      endif

!--------------------------------------------------------------------
!    calculate the flux of rainwater due to fallout.
!--------------------------------------------------------------------
      if (qrw /= 0.0) then
        dqrw3 = (qrw**1.125)*5.1*dt_micro/rr
      else
        dqrw3 = 0.  
      endif
 
!--------------------------------------------------------------------
!     calculate effect of entrainment on rain water
!--------------------------------------------------------------------
      qrw = qrw/ent

!---------------------------------------------------------------------
!    add the microphysical terms to the rain water equation. if neces-
!    sary adjust the magnitudes of the loss term so that negative
!    values of rain water are not created. define the final value of
!    qrw.
!---------------------------------------------------------------------
      qrwa = qrw + dcw2 + dqcw3 - dqrw3
      if (qrwa < 0.) then
        red   = (qrw + dcw2 + dqcw3)/dqrw3
        dqrw3 = red*dqrw3
      endif
      qrw = qrw + dcw2 + dqcw3 - dqrw3

!--------------------------------------------------------------------
!    if in diagnostic column, output the rainwater entrainment effect
!    (ent) and the updated value of rainwater. (qrw).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 2e20.12)') &
                               'in micro: ent,qrw= ',ent,qrw
      endif

!--------------------------------------------------------------------
!    apply realizability consitions to qcw and qrw. define total liquid
!    water qlw.
!--------------------------------------------------------------------
      qcw = max (qcw, 0.0)
      qrw = max (qrw, 0.0)
      qlw = qrw + qcw

!--------------------------------------------------------------------
!    if in diagnostic column, output the condensation (dcw1), the 
!    rainwater fallout (dqrw3) and the total liquid remaining (qlw).
!--------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a, 3e20.12)')   &
                      'in micro: exit micro dcw1,dqrw3,qlw= ',   &
                       dcw1,dqrw3,qlw
      endif

!--------------------------------------------------------------------
!    convert the liquid water from g / m**3 to kg / kg.
!--------------------------------------------------------------------
      qlw = 1.0E-03*qlw/rho

!-------------------------------------------------------------------




end subroutine don_cm_micro_k


!######################################################################

subroutine don_cm_freeze_liquid_k    &
         (k, diag_unit, debug_ijt, Param, tbot, ttop, qlwest, dfrac, &
          dtfr, dtupa, dfrtop, sumfrea, ermesg, error)

!----------------------------------------------------------------------
!
!----------------------------------------------------------------------

use donner_types_mod, only : donner_param_type

implicit none

!----------------------------------------------------------------------
integer,                  intent(in)    :: k, diag_unit
logical,                  intent(in)    :: debug_ijt 
type(donner_param_type),  intent(in)    :: Param
real,                     intent(in)    :: tbot, ttop, qlwest
real,                     intent(inout) :: dfrac, dtfr, dtupa, sumfrea
real,                     intent(inout) :: dfrtop
character(len=*),         intent(out)   :: ermesg
integer,                  intent(out)   :: error


      real ::  dfraca, dtupb

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!--------------------------------------------------------------------
!    define the amount of liquid in the updraft which may be frozen.
!    multiply by a factor (freeze_fraction) to take account that not all
!    this water will freeze before falling out. use Leary and Houze
!    (JAS,1980) cell_precip/cu ratio to estimate this ratio.
!--------------------------------------------------------------------
      if ((tbot >= Param%tfre) .and. (ttop <= Param%tfre) .and.    &
          (dtfr == 0.)) then
        dtfr = qlwest*Param%hlf/Param%cp_air
        dtfr = Param%freeze_fraction*dtfr
      endif

!---------------------------------------------------------------------
!    define the fraction of liquid which is to be frozen at this level. 
!    freezing is assumed to occur linearly between 258 K and 248 K 
!    (dfre), with all liquid being frozen at 248 K. the fraction is not
!    allowed to decrease from its previous value if the temperature 
!    starts to increase with height within this temperature range. the 
!    fraction of the total which is frozen at the current level is 
!    dtupb; the amount frozen at this level (dfr(k+1) is the difference
!    between this amount and the amount frozen at the previous level 
!    (dtupa). cumulative total of frozen liquid is kept in sumfrea.
!---------------------------------------------------------------------
      if (dtfr > 0.0 .and. dtfr /= dtupa) then 
        dfraca = MIN ((Param%tfre - ttop)/Param%dfre, 1.0)
        dfrac = AMAX1 (dfrac, dfraca)
        dtupb = dtfr*dfrac 
        dfrtop   = dtupb - dtupa
        sumfrea = sumfrea + dfrtop   
        dtupa = dtupb

!--------------------------------------------------------------------
!    if in diagnostics column, output the values of available liquid
!    from freezing (dtfr), the amount frozen at level k+1 (dfr) and the 
!    constants hlf    and cp_air at this level.
!--------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(4(a, e20.12),a, i4)') &
             'in cloudm: DTFR=',DTFR,' dfr=',dfrtop  ,   &
               'LATICE=',Param%hlf, 'CPAIR=',Param%cp_air,'k= ',k
        endif
      endif



end subroutine don_cm_freeze_liquid_k 


!#####################################################################



!######################################################################

subroutine don_cm_clotr_k    &
         (ntr, diag_unit, debug_ijt, Param, sou1, sou2, xe1, xe2, &
          xc1, entrain, dt_micro, xc2, ermesg, error)

!----------------------------------------------------------------------
!    subroutine clotr calculates the in-cloud tracer profiles for 
!    all of the travcers being transported by the donner convection
!    scheme.
!    author :  Leo Donner,  GFDL, 14 Jan 2000
!----------------------------------------------------------------------

use donner_types_mod, only : donner_param_type

implicit none

!----------------------------------------------------------------------
integer,                  intent(in)  ::  ntr
integer,                  intent(in)  ::  diag_unit
logical,                  intent(in)  ::  debug_ijt
type(donner_param_type),  intent(in)  ::  Param
real,   dimension(ntr),   intent(in)  ::  sou1 
real,   dimension(ntr),   intent(in)  ::  sou2 
real,   dimension(ntr),   intent(in)  ::  xe1  
real,   dimension(ntr),   intent(in)  ::  xe2  
real,   dimension(ntr),   intent(in)  ::  xc1  
real,                     intent(in)  ::  entrain, dt_micro
real,   dimension(ntr),   intent(out) ::  xc2  
character(len=*),         intent(out) ::  ermesg
integer,                  intent(out) ::  error

!----------------------------------------------------------------------
!   intent(in) variables:
!
!        sou1          in-cloud source of tracer at bottom of layer 
!                      [ kg / kg(air) / sec ]  
!        sou2          in-cloud source of tracer at top of layer 
!                      [ kg / kg(air) / sec ]  
!        xe1           environmental tracer concentration at bottom of
!                      layer [ kg / kg(air) ]
!        xe2           environmental tracer concentration at top of
!                      layer [ kg / kg(air) ]
!        xc1           in-cloud tracer concentration at bottom of
!                      layer [ kg / kg(air) ]
!        entrain       entrainment factor [ dimensionless ]
!        dt_micro      microphysics time step
!        debug_ijt     is this a diagnostics column ?
!        diag_unit     output unit number for this diagnostics column
!
!   intent(out) variables:
!
!        xc2           in-cloud tracer concentration at top of
!                      layer [ kg / kg(air) ]
!
!---------------------------------------------------------------------

!--------------------------------------------------------------------
!   local variables:

      real    ::  xeb        ! layer-average environmental value of 
                             ! tracer n  [ kg / kg(air) ]
      real    ::  seb        ! layer-average source for tracer n
                             ! [ kg / kg(air) / sec ]
      real    ::  mass_ratio ! ratio of entraining parcel mass at top
                             ! of layer to that at bottom 
                             ! [ dimensionless ]
      real    ::  delta_m    ! fractional amount of mass added to 
                             ! parcel in current model layer 
                             ! [ dimensionless ]
      integer ::  n          ! do-loop index

!---------------------------------------------------------------------
!    initialize the character string which will contain any error mes-
!    sages returned through this subroutine.
!---------------------------------------------------------------------
      ermesg = ' ' ; error = 0

!-----------------------------------------------------------------------
!    if in diagnostics column, print out an entry message.
!---------------------------------------------------------------------
      if (debug_ijt) then
        write (diag_unit, '(a)') 'in clotr: entering clotr'
      endif

!---------------------------------------------------------------------
!    define ratio of parcel mass at top of layer to that at bottom 
!    (mass_ratio) and the increase in parcel mass across the layer.
!---------------------------------------------------------------------
      mass_ratio = exp(entrain)
      delta_m = mass_ratio - 1.0

!----------------------------------------------------------------------
!    loop over the tracers transported by donner convection.
!----------------------------------------------------------------------
      do n=1,ntr             

!---------------------------------------------------------------------
!    define layer-mean environmental tracer (xeb) and the tracer source 
!    (seb).
!---------------------------------------------------------------------
        xeb = 0.5*(xe1(n) + xe2(n))
        seb = 0.5*(sou1(n) + sou2(n))

!--------------------------------------------------------------------
!    define in-cloud tracer amount at top of layer (xc2) as the sum of 
!    the value at bottom of layer (xc1) plus the amount of environ-
!    mental tracer entrained into the parcel plus the amount produced
!    by the internal source (seb). the time of insertion of the internal
!    source is given by tr_insert_time (in terms of parcel's mass 
!    increase during its traversal of the layer). if source is assumed
!    to be made available at bottom of layer (tr_insert_time = 0), then
!    seb is unmodified; if not available until completion of traversal,
!    then tr_insert_time is 1.0, and amount supplied is reduced by the
!    ratio of parcel mass between layer top and bottom. renormalize the
!    mixing ratio by the total parcel mass at the top of the layer.
!!! BUG :: seb must be multiplied by dt (perhaps as wv/dz ???)
!!!  BUG FIXED 6/2/05
!!!! BUG :: shouldn't the 1 + epm*(   ) be divided into seb  ????
!--------------------------------------------------------------------
        xc2(n) = (xc1(n) + delta_m*xeb + seb*dt_micro*   &
                              (1. + Param%tr_insert_time*delta_m))/mass_ratio 

!--------------------------------------------------------------------
!    assure that the in-cloud tracer mixing ratio is non-negative.
!--------------------------------------------------------------------
        if (xc2(n) < 0.) xc2(n) = 0.

!--------------------------------------------------------------------
!    if in diagnostics column, output the tracer mixing ratio at top
!    of layer (xc2), the layer-mean environmental tracer mixing ratio
!    (xeb) and the layer-mean tracer source (seb).
!--------------------------------------------------------------------
        if (debug_ijt) then
          write (diag_unit, '(a, e20.12)')   &
           'in clotr: xc= ',xc2(n)
          write (diag_unit, '(a, e20.12)')   &
               'in clotr: xeb= ',xeb
          write (diag_unit, '(a, e20.12)')   &
               'in clotr: seb= ',seb
        endif
      end do

!--------------------------------------------------------------------


end subroutine don_cm_clotr_k

!#######################################################################

